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1. Introduction. For turbulence produced behind the grid in a wind tunnel, Taylor’ 
introduced the assumption that the spatial pattern of turbulent motion is carried past 
a fixed point by the mean wind speed without any essential change. The assumption 
has been applied to connect the time spectrum with space correlation, and Taylor found 
indeed good agreement between theory and experiment. 

The success of Taylor’s hypothesis in the,case of wind tunnel turbulence has led 
people to use it also for shear flow. This is obviously a tentative step until justification 
can be found. In wind-tunnel turbulence, the justification of this hypothesis is essen- 
tially based upon the following two conditions: 

(a) the main flow is uniform, 

(b) the level of turbulence is low. 

One purpose of this paper is to demonstrate this justification explicitly as a basis for a 
similar discussion for shear flow, in which case neither of these conditions are as well 
satisfied. It seems plausible that an eddy of considerable size would change considerably 
by the shearing motion while it is carried past a given point by the main flow. In §6, 
estimates are made for the limitation, due to this cause, of the eddy size in the boundary 
layer, to which Taylor’s hypothesis can be applied; and it appears that the major part 
of the shear contributing fluctuations are of such low time frequencies that Taylor’s 
hypothesis cannot be immediately justified for its use in the determination of their space 
scale. 
In preparation for the analysis of Taylor’s hypothesis, it is necessary to give an 
estimate of the orders of magnitude of the various acceleration terms in the equations of 
motion. Such an estimate also gives us some insight into the nature of turbulence. The 
analysis is carried out in §§2-5. One main point is to estimate the acceleration due to 
pressure gradient. This is first done in §2 in a general manner without using Heisenberg’s 
specific hypothesis’ of the independence of the various harmonic components of velocity 
fluctuations. In §§3-4, general analyses are made of the correlation of the various accelera- 
tion terms. All the correlations are expressed explicitly in terms of double and triple 
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velocity correlations and certain two-point quadruple correlations. In §5, Heisenberg’s 
hypothesis is used in the equivalent form given by Batchelor® for correlation functions, 
and his results for pressure correlations are reproduced by a slightly different deriva- 
tion. In addition, we have calculated the magnitudes of the inertial acceleration and 
the instantaneous local acceleration. At large Reynolds numbers, they are the dominant 
terms in the equations of motion, the acceleration due to pressure forces being negligible 
in comparison. This suggests that, in a sense, the fluid particles are moving with very 
little interaction among them. Speculations based on this notion should be made 
cautiously. 

2. Pressure fluctuations. The fluctuations of pressure are related to those of velocity 
through the equations of motion. By making use of the hypothesis of independence of the 
harmonic components of the velocity fluctuations, Heisenberg obtained certain relations, 
in the case of isotropic turbulence, connecting the spectrum of the pressure fluctuations 
with the spectrum for kinetic energy. An equivalent development using correlation 
functions was later given by Batchelor. In this section, we shall develop two relations 
regarding pressure fluctuations without making specific assumptions. 

The fluctuations of velocity u; and pressure p satisfy the Navier-Stokes equations 


Ou; Ou; 1 Op 
ea = = al + v Au; , (2.1) 
ot Ox; p Ox; 
and the equation of continuity 
DU 
ox 


where p is the density of the fluid, v is the kinematic viscosity coefficient, and zx,’s are 
the coordinates of a point P. By multiplying (2.1) with the pressure fluctuation p’ at 


< 


another point p’, and taking averages, we obtain 


,. OU; l , Op ; 
(p u ) i (p aE : (2.3) 
Ox p Or; 


where we have made use of the relations 


9 ; 
(p’ sus) = 0, and (p’ Au;) = 0 (2.4) 
( 


since they are solenoidal tensors of the first rank” . From (2.3), we see that 


(2 op . aur) - = ap, —! sp.) (2.5) 
p O27; OX; \  p Oz; p O2X;/ 


\ 


(either with or without summation with respect to the index 7). We shall introduce the 


notation 
Ou; l Op 
Q=u;—, r,=--—. 
ax, p Ox, 
Then, (2.5) can be written as 
(r,a’) = (n,n), (2.6) 


which yields 
(r,a;) = (r;m;). (2.7) 
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Thus, the correlation of inertial acceleration and the acceleration due to pressure gradient 
is equal to the mean square of the latter. By Schwarz inequality, 


(r,a;)” < (wi): (ai) 
we have then 

(| w /?) < (lax |’), (2.8) 
giving an upper bound for the pressure gradient in terms of velocity correlations. It 
may be expected that 

(las |) ~ wy 


so that (| x, |*) cannot be larger than the order of ((u’))’/ d*. Thus, all the terms in 
(1.1) are limited to this order of magnitude for reasonable large Reynolds numbers, 
since, in that case, ((vA u;)*) is definitely of a smaller order of magnitude. 

For very large Reynolds numbers, one may get the order of magnitude of (| x; |’) 
by extending Kormogoroff’s concept to include pressure fluctuations. If we assume 
that 1/p? ((p — p’)”) depends, for small values of r, only on the velocity scale v = (ve)'“* 
and the length scale » = (v*/e)'’*, where « is the rate of energy dissipation, then 


; (p — p’)*) = v'F(r/n), (2.9) 


where F is a universal function of its argument. From this, it follows that 


L((2)) = 2 rv 
p Ox n 


l ae)" 7 (u*)* const. _ Vu’) 
p ((2 ~ | R, = v ; (2.10) 


This is in conformity with the result obtained above. Equation (2.10) has been previously 
obtained by Heisenberg with a specific value for the constant, but the derivation depends 
on specific assumptions. Batchelor used an equivalent set of assumptions, but he also 
pointed out that these assumptions may be questionable for small distances. Since the 
evaluation of the magnitude of the pressure gradient definitely depends on pressure 
correlations at small distances, it seems desirable that an independent argument be given. 

In the following, we shall examine all the acceleration terms in the equations of 
motion. 

3. Some remarks on isotropic tensors of rank two. We shall first develop a few 
general relations for isotropic tensors of the second rank, which will be useful for follow- 
ing developments. Robertson‘ has shown that such a tensor can always be represented 


by 


Rig = Ry(r) et + Ril) Ba , (3.1) 


r” 
where (, , & , &;) are the coordinates, and 


r=H+H+&. (3.2) 
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Thus, there are two scalar functions required for the definition of such a tensor. He has 
also shown that, for a tensor of vanishing divergence (a solenoidal tensor), there is only 
one defining scalar function. These results are generalizations of those for velocity correla- 
tions given by von Karman’ . The general form of such a tensor is 





1oF.. . raF\. aia 
Qu = Or op Site + (i + 9 2), ‘ (3.3) 


There is another kind of tensor defined by only one scalar function. This is constructed 
and perform the gradient operation twice; 


by starting with a single scalar function P(r 
thus, 
d°P(r) , 
re © (3.4) 
0&; 0 


wtr 


k 


Such a tensor shall be called a gradient tensor. By carrying out the calculations indicated 
in (3.4), we find 


1 0G : . ns aa 
Fr , eS £,£, - Goi , G = P"(r)/r. (3.5) 
Tr o7 
$y comparing (3.3) and (3.5), it is clear that a solenoidal tensor is a gradient tensor 
if both the conditions 
OG l OF . - r oF 
> == Sa G=F+-+; ai 
or oS or 2 or 
are satisfied. Then 
7 oF 
F + r— = const., 
or 
and hence 
, Y A C, ‘97 2 
F=C€, +—. (3.6) 
™ 


If we are restricted to tensors finite everywhere, the only tensors which are both 
a solenoidal tensor and a gradient tensor are constant multiples of the Kronecker delta 


[7] 


ik: 
From the number of scalar functions involved, it seems reasonable to expect that 
a general correlation tensor of rank two can be expressed as the sum of one gradient part and 
one solenoidal part. We shall now show that this can be done in a unique manner, if all 
the tensors vanish for ro. 
First, if there were two such decompositions, then 
1 ry 61) 2) ry (2) 

Rix = Pix + Qi: = Pi, + Qik ] 
and 
i — PY =" P ae 2) a. Q 1) 


tk S 


is both a gradient tensor and a solenoidal tensor, since linear combinations of such 
tensors obviously retain the original properties (cf.(3.3) and (3.5)). Thus, 4,;, can only 
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be a multiple of 6,, . However, by letting r ©, it is clear that the constant multiplier 
must be zero. Hence, the uniqueness of our decomposition is proved. 
The decomposition is effected by setting 


eee. Oe 
aed (« 5) 


(3.7) 
“Perr. 
R,=G+hH ea | 
which yields 
9 a -2 Dp 
F = ~2 41 [ a, + 3k) ar - f ™ ar, 
(3.8) 
G = —! {h [ r(R, + 3R,) dr + 2 | fi, hf 
ow J. Jr 


For convenience of reference, some of the properties of gradient and solenoidal 


tensors are summarized in the following table. 




















| Gradient tensor P;, | Solenoidal tensor Q;, 
| | 1 oF 
G = —— — 
P, =1 e468, | P= dr ar FF 

" . Yr or | 
General form 

| — | a) 

| (P’ = 1G) | + (r + 55 )ou 
Defining scalar 7 | 

P(r) =G+r oG F(r) 





(R,, for &; = (r, 0, 0)) | 








) | 
Laplacian | 2 
OF |. 40F 
—s P icity Satay die ee 
, ’ a | ’ or” + r or 
AR,, | 
. te er ar 
Trace AP = 3G+r - | 3F +r or 








If P,, and Q,, are the components of the general tensor R,, , then F and G are given 
by (3.8) and the traces of P;, and Q,, are given by 


AP = (R, + R;:) —2 [ a dr, 
(3.9) 


on Ry 
42 ('F) = +2R, + 2 | : w-{ 
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4. Correlation of the acceleration terms in the equations of motion. The equations 
of motion 
Ou Ou 1 Op 
py, t= —- Fs y ay (4.1) 
ot ” Ca p Ov; 


contain the four types of terms 


g. = ——. (4.2) 
ol 
Ou, . 
a; = Uu; > . (4.3) 
eRi 
1 Op . 
eS 4 (4.4) 
p O27; 
pw, 2 pA, (4.5) 
Among them, sixteen two-point correlations can be formed. By symmetry, only ten of 
them are distinct. Two of them, namely 
(7 a; and (7 i Me); 
are identi ally zero, § ( they can be expressed in terms of the solenoidal tensors (pa;) 


and (pu;) which are identically zero themselves. We are, therefore, left with eight in- 
dependent tensors 


By forming correlations of (4.1) at P(x;) with each of the four vectors (4.2) — (4.5) 
at the other point P’(x,), we obtain four relations among the eight tensors, thus reducing 


the number of independent tensors to four. For systematic analysis, we may choose these 
to be (a,a{), (a,ai), (r,7;) and (u,u,). The other tensors are then expressible in terms 


of these as follows 


a = (a.a;) = s{— a,ay) + (use) — Viet, (4.6) 
a.mr.) = (r.a’) = 0 1.7) 
au.) = pa} = + {(a,a; + Mh Mh — V/,,}, (4.8) 
x 7 = (rr Qty, — T af (4.9) 
an.) = wad = z{- a,al) + (uut) + Vix}, t.10) 
ri) = us) = 6 (4.11) 
where J’,, is given by 
(a a; = TT), -- F ° $.12) 


is a tensor of vanishing divergence. On the other hand 


We notice that, from (4.6), V 


(r,m;) is obviously a gradient tensor. Hence (4.12) is the decomposition of (a,a;) into 


th 


its component parts, and is, therefore, equivalent to two relations. These are given by 
(3.8) of the last section, with the correspondence given by the following relations (&; = 
, 


wy = &; 
P(r) = (pp’), 
P;, = p (xm), 
(4.13) 
Qu = p Vi ; 


R = - , 
Lik = Pp \@;Qy/. 
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Thus, (#;7r{) and V,, are both expressible in terms of (a;a;). The independent tensors 
are then (a;aj), (a,a%) and (u;uz). 

Expression in terms of velocity correlations 


We shall now show that all the above correlations can be expressed in terms of velocity 
correlations. First, let us note that, 


(ume) = v AA(uuZ), (4.14) 
(ami) = —vd k4 (ujujut), (4.15) 
0&; 
3° 
(a,a,) = 8, Ob (ujujuguy). (4.16) 


As mentioned above, (4.12) enables us to express both (2;7{) and V;, in terms of (a;at). 
Furthermore, (4.10) and (4.12) give 


(aay) = (uj;ue) — Zaspt) + (a,al) — (ry) (4.17) 


allowing (a,a;) to be expressible in terms of correlations of velocity components at a 
given instant. The other equations in (4.6) — (4.10) then give all the other correlations. 
Thus, there are altogether four scalar functions involved in expressing all the correla- 
tions of accelerations; two of them are the usual double and triple velocity correlations 
introduced by (4.14) and (4.15), and the other two relate to quadruple velocity correla- 
tions, cf. (4.16). 


Magnitude of various acceleration terms. 


By contraction and putting &; = 0, into (4.9), (4.14), (4.16) and (4.17), we may 
obtain the various quantities entering into the consideration of the magnitude of the 
acceleration terms in (3.1). We obtain 


w, |?) = (um) < (| a, |), (4.18) 
ur |?) = 350%u’) fi", (4.19) 
Op.) = +35r(u?)7ns’’, (4.20) 

' a 
a, |) = -[ x aE, (uu vu) | (4.21) 
(| ay |) = 35v7(u?)? fx” — TOv(u?)?*he”” + (| ax |?) — «| me |”), (4.22) 


where f(r”) and h(r) are the double and triple correlation functions introduced by von 
Karman and Howarth, and the subscript zero denotes evaluation at r = 0. The second 
relation in (4.18) is obtained by Schwarz’s inequality, as noted in §1. 

It is quite clear from (4.22) that for large R, = u’*/v, the mean square value (| a; |) 
of the magnitude of a, is at most of the order of {| a, |”), which may be estimated to be 
of the order of (u*)*/\*.. To make a closer analysis, one has to introduce suitable ap- 
proximations for (| a, |*) and (| x, |*). This will be carried out in the next section. It is 
convenient to put the above relations in dimensionless forms. For example, (4.22) may 
be written as 

(|e. |") ae : {| a, |*) (| we |) 
(u?)? »? _ R? (G _ R,S) + (u?)? 2 i (u?)? »”?? (4.23) 
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where 
G = fin’, S = 2hi’'r’, (4.24) 


and (4.19) may be written as 


{| ue |) _ 3564 


5S er. (4.25) 
(u A th 
5. Hypothesis for reducing the quadruple correlations. The above relations are all 
exact. To get more useful results, we shall make use of the assumption® 


uuu) = (ujupujuy) + (uyut)- (uur) + (uju; uu). (5.1) 
Then (a,;a;) can be expressed in terms of a single double correlation function. The total 
number of scalars involved would reduce to two, say f(r) and A(r). In fact, caleulations* 


show that 


(u*)"*a,at) = (4sff'’ + 14ff" — 16f”)(&e — 864) 
(5.2) 
— (A4sff’’ + l4sf” + 10ff')6,, 
where 
= df re 
= 7. f‘=—, ete. (5.3) 
. ds 
Then 
(| a, |") = 150°)" /d’. (5.4) 
By (3.9), we may obtain A P from (5.1). Thus, 
p(u’))~° AP(r) = 24 f'(8) ds + l6sf’(s), (5.5) 
which gives 
_ on ae ~« f° fat de =e 
7p |) = 24(piu’))” | f(s) ds = 12(pu')” | (4£) — (5.6) 
Jo * Jo \dy/ y 


This formula for the pressure gradient was first obtained by Batchelor. 
The magnitude of the other acceleration terms can then be easily obtained. We list 
below all these relations in dimensionless form: 


* alla 15, (0.4 
*\ eg » f (df\'d = 
. See 3(>) = 120° | ( f) <2 (5.8) 
pu r A, «0 dy Yy 


( Ou; ) 
; 30 Pe aes 2d 
eee as oe iG ~ Be + 16 — 1 | (4£) a (5.9) 


wri : dy} y~ 
er i tat 35 et ‘ 
A, |?)/{Q2)?/n?} = 120? | (‘ f) 24S eg ~- oe (5.10) 
Jo \dy/ y R) 


*'t has been found convenient to introduce s = r? as the variable in the process of calculation. 
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where A, is the total acceleration 


Ou; Ou; 
A; = at + u; az,’ (5.11) 
The symbol A, is defined by (5.8) and has been estimated by Batchelor for various 
Reynolds numbers. In particular, for very large Reynolds numbers, he found 


= 0.11R}”. (5.12) 


This is in agreement with the general discussion given in §1, and shows that the pressure 
gradient has a magnitude much smaller than the inertial acceleration. 

6. Taylor’s hypothesis. The above results will now be used to examine Taylor’s 
hypothesis. Consider first the case of a hot wire used to measure wind tunnel turbulence 
behind a grid. To the extent that the field may be regarded as homogeneous and isotropic, 
we may consider an observer moving with the mean wind speed U and reduce the problem 
to that of a registering instrument carried at a velocity —U through a field of homo- 
geneous isotropic turbulence at rest. Then, the instrument registers u(z + Ut, t), and 

i 9 

dt at Ox 
is the time rate of change as observed by the instrument. Taylor’s hypothesis that the 
instrument registers the space variation is essentially based on the argument that 


((2) ) «K uX (2) ), (6.2) 


It is seen from the above analysis that this is indeed the case provided 


To be more precise, (5.9) gives for large Reynolds numbers, 


2) (2) 5 We) 
((5i) oe) ) = 8 _ 


This may then serve as a guide for deciding on the accuracy of Taylor’s hypothesis. 
Note that the number 5 is obtained by using the specific hypothesis in §5, but it should 
be correct in its order of magnitude. 

Similar remarks apply to derivatives of higher orders. 

Another way of looking at Taylor’s hypothesis is to say that the main terms in the 
equations of motion 


me 41 u Se +, in a 2 +» Au, (6.4) 
are the first two, the others being much smaller in comparison. This is seen from the 
rough argument that the pressure terms are of the same order or smaller than the con- 
vective terms u,du,/dz; as justified by the previous analysis. Indeed, the pressure terms 
are often much smaller than indicated by such a crude argument. 

For the case of shear flow, it is difficult to make an exact analysis, and one may first 
proceed with the kind of rough arguments just given, with its plausibility supported by 
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its correctness in the case of uniform flow. However, the pressure fluctuations are not 
only related to the quadratic terms in the fluctuations but also to the convection of 
momentum of the mean motion due to turbulent fluctuations. Thus, if U(y) is the basic 
parallel flow, in the equations of motion 





Ou , Ou dU Ou 1 Op 

a l = = ss oe = ot + vAu, 

Ot Ox dy Ox; p Ox | 

Ov , Ov Ov 1 Op — 
— + U- +u,— = —-£4+ vav, | (6.5) 
ot Ox " OX; p oy 

Ow , OW ow l dp 

—— | [] — hh = = oF + vAw, 

Ot Ox ' OF; p Oz 


there is the term v dU /dy. If Taylor’s hypothesis were true. Then u = f(x — Ut), v = 
g(x — Ut), w = h(x — Ut), and we have approximately, for large R, and small turbu- 


lence levels, 


ap __ aU | 
Ox - “a dy’ 
Op ; —r 
=e, | (6.6) 
OY | 
| 
9 | 
— = 0. 
Oz ) 
Hence, 
1U , 
a= —% s Ga - Ud, where G’(é) = g(é), 


and the second equation of (6.6) gives 


Op ro 4. (uy ~~ gf =— +t 
= ou» eee G - — « — +} ~ -4 
oy p dy ‘s) p dy g's) 0 


this is generally true only when 


dU _ 0, qU =0 
dy 


dy — 
i.e., when the main flow is uniform. There is, therefore, no general justification of extend- 
ing Taylor’s hypothesis to the case of shear flow. 
One may still attempt to justify it on the basis that 
ou dl 


sl la 


Then it may be expected to hold only for components of wave length k such that 


Uy Wu ay yy ¥ WY 
kU > a or ky > U dy’ 
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The latter form is useful when the experimental data are presented in scales of log y, 
For boundary layers, experimental values at y/6 = 0.2 (Townsend’*) give roughly 


ky > 0.9 or ké > 4.5. 
In terms of the wave-length A of the disturbance, this means that 


A “a 
=e Cs. 
6 
Thus, the wave-length must be much less than the boundary layer thickness in order 
that Taylor’s hypothesis may apply. 
For smaller values of y/6, the restriction is even more severe. For larger values 


y/6, the condition 


ky> 1 


may be taken as a rough approximation, until one comes to the region where inter- 
mittency becomes appreciable ( = 0.02 in Townsend’s experiment). Before this point 
is reached, and with y > 0.2 6, y(dU/dy) and U both increase slightly with y, making 
the resultant change in the quantity (y/U) (dU/dy) fairly small. 

Thus, at y/6 = 0.6, we must have 


ké > 1.5, 
Ol 


A 
=< * 
5 <K 
This means that for harmonic components with wave-lengths of the order of magnitude 
of the boundary layer thickness, Taylor’s assumption is barely applicable at y/6 = 0.6. 
Klebanoff and Diehl’ presented their results in terms of nL,/U, where n is the time 
frequency in cycles, and L, = 0.4 6. The above condition becomes then 


nhs 0.1. 
l 
With reference to their Figure 25, it is seen that this condition is satisfied from the 
n-*’* range upwards. Shear-contributing fluctuations, however, do not satisfy this 
requirement. 

It should be noted, however, that the wave-length U/n that may be assigned to the 
n-** range is of the order of magnitude of the thickness of the boundary layer (6/2 
to 6). This is rather large for the application of Kolmogoroff’s theory. The possibility is 
not to be excluded that the disturbances are actually propagating downstream at a 
lower speed analogous to oscillations in a laminar boundary layer. 

Referring to (6.3), it is seen that for the observed level of turbulence in the boundary 
layer, the non-linear terms does not cause much inaccuracy in the application of Taylor’s 
hypothesis. This is contrary to the conjecture that might be based on the fact that the 
quadratic terms are important in the phenomena of turbulent motion. 
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NON-UNIFORM SUPERSONIC FLOW* 


By 
A. ROBINSON 
University of Toronto 


1. Introduction. The forces which act on an aerofoil in a non-uniform main stream 
have formed the subject of a number of papers (refs. i-5). In all these investigations, 
incompressible flow only is considered. References 1-4 are concerned with the effect 
of non-uniformity of a two-dimensional main stream, while reference 5 deals with the 
effect of non-uniformity in spanwise direction. Much of the work mentioned was done 
in connection with wind tunnel problems, but that application may be less important 
in supersonic flew. However there are other cases of non-uniform supersonic flow which 
may be of practical interest, such as the flow around jet vanes.’ In any case, the physical 
significance of the problem appears to warrant its investigation. 

In the present paper, we shall analyse the problem of accelerated supersonic flow 
for the two-dimensional case on the basis of linearised theory. According to that theory, 
the effect of any small curvature in the main flow would be regarded as equivalent to a 
curvature of the aerofoil in opposite direction and there is therefore no need to consider 
the problem further. To be sure, acceleration will be associated with curvature some- 
where in the main stream, but in order to consider the problem of acceleration in isola- 
tion one may suppose that the aerofoil has been placed in a plane of symmetry of the 
main flow so that the velocity component of the main flow in a direction normal to 
that plane vanishes throughout. 

The acceleration of the main flow can be measured conveniently in terms of the 
variation of the static pressure. We shall derive closed analytical expressions for the 
pressure distribution around the aerofoil in the particular case of linear variation of 
the static pressure in the main stream. Simplified formulae will be given for sufficiently 
high Mach numbers. These do not apply to Mach numbers near unity, but in that 
region the entire theory, like the theory of uniform flow, will be less reliable in any case. 

2. Linearization. As stated in the introduction, we assume that the aerofoil is 
placed, approximately, in a horizontal plane of symmetry of the flow, and we take the 
z-axis along the direction of the main stream in that plane, and the y-axis in a direction 
normal to it. Let U(2,y), V(x,y) be the velocity components of the main stream, and 
U (x,y) + u(ax,y), V(2,y) + v(2,y) the velocity components in the presence of the aerofoil. 
Both the main flow and the flow induced by the aerofoil are supposed to be irrotational 
(at least within the approximation adopted) so that 


OV aU 

. —— =o ( 2.1 

Ox oy . 2.4) 
and 

Dn Fn (2.2) 

Ox oy 


*Received Nov. 23, 1951. 
1This remark is due to a referee to whom I am indebted also for some further suggestions. 
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We infer the existence of an induced velocity potential ¢(z,y) such that 


od 0 oi 
“ua ==, y= = J (2.3) 

OX oy 
In addition, the equation of compressible flow must be satisfied both in the presence, 


and in the absence, of the aerofoil so that 


°\ aU V*\aV 2UV aU 
(, - Ui) aU , (Hav _ aur aw g - 
as Ox a/ oy a oy 
and 
(U +u)’\ aU +0 ( (V + a) av +) 
(1 a” Ox ry a OY 
(2.5) 
_,U+twV+y)aU+y _y 
= a” oy 7 
where a” = dp/dp is the square of the velocity of sound, p being the pressure and p the 


density. Strictly speaking, a is a function of (U + u)? + (V + v), but we shall assume 
that the proportionate changes of pressure and density are small so that a may be re- 


garded as constant. 
Expanding (2.5) and taking into account (2.4) we obtain 


((1 — &) - (2H + 4)) U — (Pe 4H) ey (1 - BY 
a’/ \ as} dx a a/ 0x a’/ oy 


a 


2Vv % ar fi QUV du 2Uv . 2uV ap) au ou) 
= ae a a — a 53 a ee ee —— em: lead a ar i sz = 0 
( a® a/\ oy OY, a oy ( a” + a + a /\oy + Oy 


(2.6) 


We suppose that U is of the same order as a, so that U/a is of the order unity, but 
that U/a differs from unity numerically to the extent that 1 — U’/a’ also is of order 
unity. Furthermore, we suppose that V, uw, v are small compared with U and a. Regarding 
(2.6) as a linear form of the derivatives du/dz, --- , dV/dy and omitting from the co- 
efficients terms which are small compared with other terms still included, we obtain 

( _ 2) du _ 2UuaU | ww _ 2¥v+v°aV _ 2UV du _ 2Uy (3z au) _ 
a’/ dx a” ox * oy a” oy a oy a \dy ~ oy 

Moreover, although that assumption may break down locally, we suppose that the 

relative orders of magnitude laid down above persist for the derivatives, so that we 


may omit the terms 


2Vu+eaV = =2UV du q 2Uu du 
a° oy’ a dy’ . a” oy 
Then 
(1 — 2) — Seed , we _ Bea og (2.7) 
a/ Ox a ox oy a oy 


Finally, we suppose that the variation of U in the direction of the y-axis is small 
compared with the variation of U along the chord, dU/dy «K dU/dz so that we may 
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regard U as a function of z only. Since 9U/dy = dV/dx = O in the plane of symmetry, 
this condition must be satisfied at least in the neighbourhood of the aerofoil. In con- 
sequence, we shall neglect (2Uv/a’)dU/dy while retaining (2Uu/a’)dU/dx. Equation (2.7) 
then becomes 





vu) du 2UdU ov 
(: a) ax a dx” A dy . 


Or, in terms of the induced velocity potential, ¢, 


Ue 1)% 2U dU a6 _ a* _ 
( 3 dx” ' a? dx dx Oy 0. (2.8) 


Let a = a(x) be the local incidence at a point on the surface of the aerofoil, taken 
as positive when the slope dy/dz is positive. Then the boundary condition is 


V +0 = (U + w tan a(2). 


We shall assume, as usual, that this condition is satisfied at the projection of the 
point on the z-axis so that V = 0. Hence, for small a, 


v = U(az)a(z) (2.9) 


at the aerofoil. 
To calculate the pressure, we have Bernoulli’s equation, 


[P+swt+or+ ty) =H (2.10) 


where H is constant throughout the medium. Or, neglecting terms of second order of 


smallness, 


[P+ 5 + 20m =H. (2.11) 


Let po , po , Uo , Uo be pressure, density, and longitudinal velocity components at 
an arbitrary but fixed point. The proportionate changes of p and p are small, by as- 
sumption, so that if 

p = pl +8), 
then p = kp” = kpo (1 + ys) = po(1 + vs), dp = poy ds, and so 


[se ee fas, pe _ lt) p= 
Po Jo L+s ‘ 


‘pe P 





Po Po Po 


Accordingly, (2.11) may be rewritten in the form 
p — Po = —5 [(U* — Ud) + AUu — Usw)). (2.12) 


It should be observed that we cannot now neglect Uu compared with U” since the 
difference Uu — U up is not in general small compared with U? — U3. 
At a point upstream of the wing u, = 0, and so 


P— Pp = 3 (U? — U2 + 2Uu). (2.13) 
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When the aerofoil is absent, 


2 


Dp — po = ——[U? — Us, 


so that the static pressure is a function of x only, within the approximation adopted 
here. More particularly, we shall assume that the main stream static pressure is a linear 
function of z. Assuming that po , po , Up correspond to the origin, we then have 
Dp = Po + Qe. 
Hence, using (2.13), 
9 9, 9 
. - 2 2 2g ey 2 
U=U2?+=(p—p) = U2 -Sx=U24+r, d= ——S 
p p Po 


Thus, \ > 0, g < 0 for accelerated flow and \ < 0, qg > O for decelerated flow. 


3. Accelerated flow. Let JJ = U/a be the Mach number of the flow, and 6 = 
V M° — 1, so that J/ and @ are functions of x only. Then (2.8) may be rewritten as 
97 ] a 4 4 ae , 
SE + = (6?) SSF = 0, (3.1) 
Ox dx Ox Oy 


Taking first the case of accelerated flow, \ > 0, we introduce a new variable r by 


2a? 2a’ ) > DN 
r= -f = = Bo + 32, » > OD. (3.2) 

A A \ \ ee -- ; 

Then 
Ar’ a” 3? 
.= 5 Bos 
4a” A 

and 

dr | ; 2a" 

dx hh mM 

\ Bo +t a4 
a 
Hence 
dd _ dA dr _ 2a" 0 
0x0 Or dx”—s wr: Or’ 
Also 
d d (Us + rx r 
fm (a - 1) = 5, 
dz dx \ a a 

and so 


dp _ 2a° a ( 2a" ae _ 4a‘ a> _ 4a° dg I (zs _ 1 26) 
On" Ar Or \dr Or] *r’ Or” Nr? or B° \ar" r or/” 


J 


Substituting in (3.1) 


9 10 a°¢ sr 
ce ht ae a Soe eG, (3.3) 
Or YT or oY 
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Now consider the partial differential equation 


20a (3.4) 


This may be regarded (amongst other things) as the linearised equation of steady 
supersonic flow for a Mach number ~/2, where the main stream is directed along the 
f-axis (compare refs. 6, 7, 8). Particular solutions of (3.4) are obtained by means of 
source distributions over the ¢, ¢-plane, — > 0, thus, 


¢ a a o(é, ’ fo) dé, Ago —— ~ 
o(E, Y, 9) ~ Se I Vi — &)* ae y oe iF = bo)? (3.5) 











where the area of integration on the right hand side is given by (§ — &)? — y? —(¢ — 
tf)” > 0, & — & > 0. The source strength o(f, f) is connectéd with the normal deriva- 
tive of r at the é, ¢-plane by the relation 


rs) 1 
oo a a G6) 
In (3.4), introduce new variables r, Y by means of 
§ =rcosh y, 7 =rsinh y. (3.7) 
(3.4) then becomes 
do , 106 _ 10% _ d% _ 
or ' rar roy oy | % (3.8) 


while (3.5) may be written as 


¢(r cosh y, y, rsinh y) 
1 I a(ro cosh Yo , To Sinh Yo)ro dro dy 
Qe Jd 





Vie cosh ¥ — ro cosh ya)? — y? — (rsinh ¥ —resinh YF, 

or 
¢(r cosh y¥, y, rsinh y) = 5 I = 8 PT. Soe Se Me ; (3.9) 
2r Vr +7 — 2rr cosh (¥ — Yo) — ¥” 








where the area of integration is given by 
r cosh ¥ — ro cosh Yo > 0, r+rn- 2rr, cosh (¥ — Yo) — y > 0. 


Suppose now that ¢ is independent of Yo , o = f(r). 
Then (3.9) becomes 
¢(r cosh y, y, rsinh y) 
= x | fran ar | avo — 
an eS Se + 22 = Dre cosh (¥ — yo) — y’ 


We denote the integral with respect to y on the right hand side by J. Substituting 
w = ¥ — wv in that integral, we obtain 





dw 


jaf , —_ 3.11 
Vr? + 735 — 2rro cosh w — y° (3.11) 
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where the range of integration is given by 


r? +r — 2rro cosh w — y? > 0 
wi 7. Nast (3.12) 


Hence J is independent of ¥, J = J(r,y,r0). If follows that the function ¢, which is given 
by (3.10), i. e., by 


¢= os [ J(r, y, To)f(To)ro aro (3.13) 
™ J 


is independent of ¥. Nowsall integrals (3.9) satisfy (3.8), and so ¢ as given by (3.13) isa 
solution of (3.3). Moreover, by (3.6), 

26) ” 

— = —— {(7,). 3.14 

(32 ai ail y) (To, ( ) 

This suggests that we may solve (2.8), for the case of uniformly accelerated flow 

and for the boundary condition (2.9), by a suitable choice of f(79). However, in addition 
to satisfying (2.9), the physically correct solution of (2.8) must vanish ahead of the 
Mach lines (characteristic curves) which emanate from the leading edge of the aerofoil 
in a down-stream direction. The equation of the characteristic curves of (2.8) is 


(Z _ 1) dy’ — dx? =0 (3.15) 


so that the characteristic curves which emanate from the leading edge in downstream 
direction are given by 

dy U’ ‘ise 

oY (Z _ 1) for y > 0, 


dx F 
(3.16) 


dy _ (SZ y" : 
fi a 1 for y < 0. 
The corresponding curves in the r, y-plane as given by (3.2) are characteristic curves 
of (3.3). They have the equations 
y = +r-+ const. (3.17) 
In particular, if the leading edge of the aerofoil is at the origin of coordinates in the 
xz, y-plane, then the corresponding value of r is r = (2a°/A) and so the two Mach lines 
in question are given by 


2a’ 

y = +(r — 7’) where = “7 Bo - (3.18) 

To satisfy the condition that ¢ vanishes upstream of these curves we only have to 
assume f(r,) = 0 for 7 < 71’ or 


: 
¢ = = | J(r, y; Tf (Toro Aro « (3.19) 


r? 
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To carry out the integration in (3.11) we introduce the variable 





2 dt 
t=cosh=, dw =—=———. 3.20) 
2 ” Vt? ag ( 
Then 
r+ ri — 2rr, coshw — y’ = (r +7)? — y’? — 4rro cosh? 5 
= (r+) - v(t — k’ cosh? “) 
where 


_ __— Arr 


= (r rn ro)" aa y*" (3.21) 





We observe that, for r > 0, 75 > 0, 
(r + 7%)? > 4rro 


so that, given 7, 7 , the right hand side of (3.21) is positive but not greater than 1 for 


sufficiently small | y |, and we may assume 0 < k < 1. For such y, (3.11) becomes 





4 7 [~ dt 





J = SS oem 
/r+ner—-—yn Vl — DO — ke’) 
(3.22) 
— . - K'(k) = —— s some, K(k’) 
Vir+n) -—y Vir+nrn) -y 


where K, K’ are the complete elliptic integrals of the first kind, and the complementary 
modulus, k’, is given by 
41m) (r =f)" oe y 


 r+tn)?— y (r+n)-y 


o 


kb” = I 
In particular, for y = 0, 
4 jr? 4 41 — (f/r) 
r= Aas) staal ste) 
r+ Ay +1} T+ AMY + (r0/r) 
But 
K(L= ro/t)) ~ (1 4 ro) ("2) _t+% x’(") 
‘\T + (r/r) ~ r r})~  & ™\y 
by Landen’s transformation, so that, for y = 0, 


J= 2 K("), (3.23) 


ry 
Substituting this expression in (3.19), we obtain 


f(r.) 22 x’("2) dr, (3.24) 
Tr rT 
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The function f(r) is given by (2.7) and (3.14), 


i) = 20a. (3.25) 
Now 
os i 2a — 
de a 1 x VM* -1 by (3.2) 
and, for any modulus k, K’(k) = K (W1 — k°), so that (3.24) may be replaced by 
4c’ ¢* a { [M? —m\ , 
ie ae (ST) am 

where 


(3.26) 
Me = 5/4. M = 
4. Approximate treatment. 


U(x) /a. 


The elliptic integral K in (3.26) can be expressed in 
terms of the hypergeometric function F, as follows 


: Mw m T (2 ] M? — m’? 
K(4 M’ —1 ) ey i + | 
so that 


99°) PF — 4 ) (4.1) 
a(R) = 5(i+ LR) am 
except for higher powers of (J? — m’)/(M’* — 1). With this approximation 
¢ = 2a. wal = [" ~ m? dm + avi <n [ a(m*? — m’*) am} (4.3) 
and so 
dd do dM 
Oe an SOOM a 
on meal as ae oe ‘wart i ; am dm 
VM -1 VM? —1 4 — 1) Jy. 
3 a 


_- - > oe a( , — 1 } (4.< 
ir — D? ‘ m m) dm 4) 


Assume now that a is constant in the interval 0 < x < c¢. For points in that interval, 
the expression in the curly brackets on the right hand side of (4.4) equals, approximately, 
4 2 9 ( Po 4 
(_SeM*_ , Sota — 1) 

4(M*—1) °° 4(M* -— 1° 


1) 


Hence 


. 1 M? 
(M — M,) = 9% Wk A (M —- M,). 


(4.5) 
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where 


_ M(M — Mo) 
h= mM? — 1) ° (4.6) 


We observe that — U(x) a/ VM’ — 1 would be the value of wu if the theory of uniform 
flow applied at the point under consideration, so that h may be regarded as a correction 
term. The approximate formula (4.5) applies only provided MZ — M, is small compared 
with M and M, , and provided these Mach numbers are not too close to 1. For example, 
if M, = 14, M = 1.6, then h = 0.10. 

An alternative expression for u is obtained by writing (4.5) in the form 


M M’ ) 
u= aa( 4 - XM? — 1)” (M — M,)}, 





and expanding the right hand side at M, . Then 


M, 2 + M; ) : 
= “at. , 6a ye ~ 
l aol VMs -— 1 2(M; — 1) ( 1 1 (4.7) 





except for higher powers of M — M,. 
The following expressions for the pressure are then obtained by substituting (4.5) 
and (4.7), respectively, in (2.13). 


Ua 


p = po — 2 (U? — Us) + po Sara = (4.8) 
and 
Po 2 +2 Usa r M; — « 
= - 372 _— : - ——neng ners (Mf — M,). 9 
p I 9 6b Us) + p Vz a + PU oe oc? — 1) VM,) (4.9) 


It will be seen from (4.8) and (4.9) that for given constant a, p depends only on the 
local Mach number and on the leading edge Mach number and (free stream) pressure 
but not on the absolute rate of acceleration. In other words, p is the same for points z 
and x’ = u» x on two aerofoils at incidence a if the main stream velocity and free stream 
pressure at the leading edge are the same in both cases, while the rate of change of the 
square of the main stream velocity, d/dx(U*), along the second aerofoil equals 1/ times 
the rate of change of the square of the main stream velocity along the first aerofoil. 
The same conclusion is reached if we express the pressure in terms of the more exact 
formula (3.26) for the velocity potential. 

Assume now that the local incidence at the top surface is a for 0 < x < c/2 and —a 
for c/2 < «x <c. For0 < x < c/2 the pressure is then still given by (4.8) or (4.9), while 
for c/2 < x < c we now have instead of (4.5), 


Zax 


one = 
= pq) (4.10) 


| 


Uu 


— 


where 


yr. MIM — M,) — (M, — M,)] 
ai 2(M* — 1) , 








(4.11) 
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M, being the Mach number at x = c/2, M = U(c/2)/a. 
Similarly, (4.7) is replaced by 


s 





M, M2 +2 ; M? | 
“= aa(— eS —— 5 a. a ee (M — M.) + aa pa A M,)). (4.12) 
a. @\stL oO 4m£ 0 


Accordingly, the two expressions for the pressure are now 





Po 72 72 U*a ‘ 
paz kh - (U ver l o) — po” 7s — (1 — h’) (4.13) 
pee 2 VM* -1 
and 
er ; Ura - M; — 4 
=) — f ‘Vek i eo = a os ae a = . 
Pp Pp 9 ( 0 Po V/ M es 1 poll o@ 2M (MM? — 1)” (M M.,) 
(4.14) 
= M, : 
— poUoa OW? — ))* (M, = M,). 
5. Decelerated flow. For decelerated flow, \ < 0, we put 
2a° sf f. . k | 
r= —-—f=—-—4/+az, r>o. (5. 
, 4 x Yh + 22 r>* (5.1) 
Then 
_ A » = a’ 3? dr - 2a" 
~ 4a? 5? dx» 


as before. Thus the differential equation (3.3) still governs the flow, and we may adopt 
the same procedure as in section 3. However, r now decreases with increasing x and so 
the characteristic lines in the plane which correspond to the Mach lines (3.16) in the 


physical plane are 


=r —7 for y > 0, 
5.2) 
y= —(r’ — r) for y < U, 
where 7” = — 2a°/d ® . Accordingly, we now have instead of (3.19), 
<= J(r, Ys To)f (To) To Aro (5.3) 
aa Je 


where J(r,y,7o) is still defined by (3.11). Carrying out the integration in (3.11), we have 
to take into account that r is now smaller than or equal to 7, . We then obtain, for y = 0, 


2 ut 
J = —  K’ (5.4) 
ri To 


and so 
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Using the approximation 


( [m?—M x 7 _ : 
Moar) $0439). - 


and adopting the same procedure as in section 3, we find the same approximate formulae 
for the longitudinal induced velocity and the pressure as given in section 3 for accelerated 
flow. 

6. Lift and drag. Consider now a thin flat aerofoil of chord length ¢ at incidence 
(Note that according to our definition, a is positive when dy/dz is positive.) Then the 
pressure is given by (4.8) or (4.9) for points on the top surface and, by familiar considera- 
tions of symmetry, by 


Po ;y72 72 U*a 
= “Bae l a Us = oO me (t = h 6.1 
p=p . ( i, Vi oi 0) (6.1) 
or 
" —— ‘ty a Tee tok. 
DP = po — pol Us) — Po 5, eam poUoa 2M,(Mz — 1)°*” (M — M,) (6.2) 


for points on the bottom surface of the aerofoil. Hence the pressure difference between 
top and bottom surfaces is 











_ 2pUa _ ‘ 
Ap = =e (1 — &), (6.3) 
i ret + 
or 
U ca r2 _Ms smn 
= 2p) —— vo .+ poVoa Wray (M2 — 1)” (M — M,). (6.4) 
Using (6.4), we obtain for the resultant lift 
Be 1 4a . Mo-—4 _ 
L= as | Ap dz= ~9 poe oc - JM — a po Uoa M (M2 mae 1)” 2 [ (M M,) dx 
4a M; aR 4 2M, + M, (M, pe M,) 





] _ en 
= 9 oer /M? — 1 pola M.(M2 — 1)°*” 3(M, + M,) 


where M, is the free stream Mach number at the trailing edge. Replacing (2M, + M,)/ 
3(M, + M,) in the last term by 1/2, in harmony with the approximations adopted so 


far, we obtain 


4a M, — 4 ) P 
_ ———_—_—— ———,——__ (M, — M,) }. f 
L 5 poU3 of /Mm — =(1 + 4M,(M, ai y ( I, 1.) (6 D) 
Thus 
4a 
C,= al + n) (6.6) 


V Mo — 
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where 


Mo — 4 ‘a 
n = ——~,— (Mz — M,). (6.7) 
1” 4M (M3 — 1)? 
The correction term, 7, changes sign for M, = 2, and it is easy to show that for 


larger values of Mo, JJ) > 2, 7 remains negligible for reasonably small values of IJ, — My, 
| M, — M,| < 0.2, say. 
The wave drag corresponding to L is 


D = —La. (6.8) 


Next, we assume that the top surface of the aerofoil is at incidence a for 0 < x < ¢/2, 
} 


and at incidence —a for c/2 < x < ¢, while the bottom surface is at incidence —a for 
0 < xz <c/2and at incidence a for c/2 < x < c. This is the case of a symmetrical double 
wedge (diamond-shaped) aerofoil. The pressure at both top and bottom surfaces is then 
given by (4.8) or (4.9) and by (4.13) or (4.14) for the front and rear halves of the aerofoil 
respectively. The wave drag of the aerofoil therefore is 


D= | | pa dx — | pa a) 


M,) 


Il 
.~) 
~ 
| 
~ 
i 
+ 
~ 
— 
l 
+ 
a | 
—) 
bee 
— 
‘Sl hen 
|S 0 
~ 
te 
| 
| + 
_— 
= 
~~ 
ee 


re" a 


9» @Q M, 
+ poU'ce ~ 


2 (Me — 12 As — Mo) 


where U is the free stream velocity at the trailing edge, U = M/, a. Now 


M,-M,.M?—-M 1 


M.—-M,' M@-M~ 2° 


i \ 4a? 3M; -—8 ) , 
= Pa pe l1]}/+---— i . ] : (6.6 
om) E \Mi *) © Vt — 4 (+ aap Me ~ Mo sine 


and so 


\ 


a (Mo _ 1) — 3Mo — 8 


. = Tr (M, — Ma)). (6.10 
2\M Mtn i \ * aaa —p*~ * 6.10) 


We observe that in this formula a may also be interpreted as the maximum thickness- 
chord-ratio of the aerofoil. More precisely, the latter is r = tan a. 
For a positive angle of attack the lift on this aerofoil is the same as for a flat aerofoil, 


while the wave drag is the sum of the drag at zero incidence as given by (6.9) and of the 


f 


drag associated with the incidence (6.8). 
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A PROBLEM OF FINITE BENDING OF TOROIDAL SHELLS* 


BY 
R. A. CLARK AND E. REISSNER 
Case Institute of Technology Massachusetts Institute of Technology 


1. Introduction. In the present paper we consider a non-linear problem involving 
finite axi-symmetrical deflections of toroidal shells with circular cross section. The 
problem is that of an expansion joint for two straight sections of a cylindrical shell 
loaded in axial direction (Fig. 1). The particular kind of expansion joint investigated 
is generally called an “Omega” joint. 
































= an anal l -< a r r 
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Fic. 1. 


The linear theory of small deflections of such joints was recently considered by the 
first-named author [1]. The purpose of the present work is to establish the range of 
validity of the linearised theory and to obtain appropriate corrections to the linearised 
theory if the behaviour of the joint is slightly non-linear. 

The results presented in what follows may, of course, be extended with relatively 
little modifications to other loading conditions for toroidal shells. The basic differential 
equations for finite axi-symmetrical bending of thin shells of revolution which are used 
in what follows may be found in a recent paper of the second-named author [2]. 

For the present problem, under assumptions which are discussed further on, we have 
to deal with the system of non-linear differential equations 


8” + uysint = wQ(cosé + Bsin £) + wy6 sin &, 


I 


—}up" cos &. 


y’’ — pBsiné 
*Received January 2, 1952. The present paper is a report on work done under the sponsorship of 
the Office of Naval Research under Contract N5-ori-07834 with Massachusetts Institute of Technology. 
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subject to the boundary conditions B(+2/2) = ¥(+72/2) = 0. The constant parameters 
are u = [12(1 — »’)]'’(b?/ah) where b is the radius of the cross section of the middle 
surface of the shell, a is the radius of the center line of the torus, h is the wall thickness 
of the shell and » is Poisson’s ratio; and 2 = [12(1 — v*)]'’?(aP/Eh’) where P is the axial 
load per unit of circumferential length and E is the modulus of elasticity of the shell. 

As for the linearised problem, it is found that the following two ranges should be 
considered separately: 

G) #2 = OC): (iil) p> 1. 
When yz = O (1) the linearised theory is applicable as long as pQ < 1. 

When yu > 1, which is the case for many practical applications, the linearised theory 
is applicable as long as p*’* 2 < 1. 

In these two cases non-linear corrections to the results of the linearised theory may 
be obtained by developments in powers of »Q and yu“, respectively. The present paper 
contains numerical results for the second, less simple, of the two cases. 

The physical significance of these results may be seen by considering the relative 
ris deflection 6 of the two edges of the toroidal joint. As long as » = O (1) we have 

= O(buQ) = O(Pb*/Eh*) and the condition for linearity is 6 « b. When pp > 1 we 
si § = O(b2) = O(Pba/Eh’) and the condition for linearity may be written in the 
form 6 < b/z’”*. In the former case the relevant order of magnitude conditions are the 
same as those which hold if instead of a toroidal shell with center line radius a we had 
a straight cylindrical shell with a = o. In the latter case the finite radius a affects the 
results and non-linearity must be considered for progressively smaller values of 6 as yu 
becomes larger. In either case the deflection 6 may be large compared to the wall thick- 
ness fA of the shell without the occurrence of non-linear effects. 

2. Differential equations for finite axi-symmetrical bending of thin toroidal shells. 
Let 

Tr =a-+t bsiné, 2 = —bcosé (1) 
be the parametric equations of the undeformed middle surface of a toroidal shell of 
constant thickness h. The basic differential equations for finite bending, up to terms 
of second degree, are then [2; Eqs. III and IV], 


, 4\2 ? € A vr 
Te r6 To 3 Talo Vv Zo ° 
Ar am , “a ae i : a ool 
p’ + 2p - (8) -. p-|5 72-22 |e 
To To To a 15 2To 


<7 {W sin — (ro V) cost — B[V cosé + (7, V) sin £]}, 


py 4 ry! | (#) 4 yo hk us [2 28 , , 2 ay _ » = g ay 


b c; ° 9 pt 4 2, 
= — E sin § — 58° cos | + E + = a Vet y+ ‘ TV)’ (3) 


0 0 


(zi)? — (ri)? ro’ |, ro 7) 
+ (Sea _, BroV) — »™ [8(roV)V 
To To 0 


br, + bz6B,.2 . 
~ (ropu)- 


(2) 


t 9 
. - (ropu)’ — 


0 
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In equations (2) and (3), V is a stress function, 8 the change in the meridian angle due 
to deformation, p, the horizontal component of surface load intensity and V the vertical 


(axial) stress resultant given by 
foV = -| (a + bsin £)bpy dé, (4) 


where py is the vertical component of surface load in intensity (Fig. 2a). The quantities 
D and C are defined by 


Eh’ 
= = C = Eh. 5 
Ds oy Eh (5) 


Primes indicate differentiation with respect to &. 
In terms of ¥ and 6 we have, for the horizontal (radial) stress resultant H 


H=WV/n, (6) 

and for meridional and circumferential stress resultant N; and N, , as indicated in Fig. 2b, 
roN; = WV cosé + (r,V) sin € + B[W sin —€ — (ro V) cos &], (7) 

bNg = WV’ + bropn .« (8) 


The meridional and circumferential stress couples J; and M, are given by 


bM, = pa (2 3 + :2 ) | (9) 
bM, = p| (Ha + 32 at) ty 2’. (10) 
ri Vv 
N 
5 








Fig. 2a. 
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Fia. 2b. 


The horizontal and vertical components of displacements u and w are given by 


u = (N, — wN), (11) 


reer 
w= —b [ (s cos — + 5 8 sin r) dé. (12) 


Omission of all non-linear terms reduces equations (2) to (12) to the corresponding 
equations of the linearized theory which were considered in [1]. 

3. Dimensionless form of the differential equations for a shell subject to edge loads 
only. Let P be the magnitude of the axial edge load per unit of circumferential length 


measured along the section = —7/2 of the toroidal shell and assume that no distributed 


surface loads are acting. This means that in equations (2) to (4) 





Pu => 0, Po = 0 (13) 
and 
rV = aV, = (a — b)P = const. (14) 
As in the linearized theory one may define two parameters \ and yp by 
i = D/a, mw = [12(1 — v*)]'7(b’/ah), (15) 
and introduce a dimensionless stress function ¥ and a dimensionless load function Q in 
the form 
(12(1 — »’*)]'” 2i—r)]'? . 
Ss Q = —-—;—“*—- eas 16 
Ele Ele ' (16) 


With this the differential equations (2) and (3) assume the form 


4 C08E oy _ a) vd sin £ | 
yt 1+ \sint? (; + rsiné + 1+Asiné B 


2 . 
Aa 3 A sin cos & _ l : vr cos & 


2(1 + Asin £é)’ 1 heat | ug) 


= ——_+ [y(sin & — B cost) — Q(cosé + Bsin £)j, 


1+Asinég 
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or. | heed)? _ shin | 
v Ts ieat’ (ae 1+d\anel” 


~ | 2 Asin boone, 5 vr cos — |v - vd sin — a’y 








(1+ sing)’ 1+Asiné 1+Asin¢ 
a ae ee 
“4 tang | esine 5 8 cos (18) 





4 of d’ sin £ cos — vd cos £ 
(1+ Asiné)* © 1+Asiné 


A*(sin’ £ — cos" é) vAsing | vr cos~ ., 
ine 1+ Asin=? P 


(1 + Asin é)’ 1+ Asiné 

We note that the definitions (15) for \ and yw and the fact that b/h > 1 for a thin 

shell always imply that \ < yu. In what follows we restrict our attention to problems 
where we also have 

\<1, (19) 

which means that the radius of the cross section of the torus is small compared with 


the radius of the center line of the torus 
Under the assumption that (19) holds the differential equations (17) and (18) can 


be simplified considerably so as to read 
8’ + pwsinéy = u[Q(cos~— + Asin £) + ¥6 cos £], (20) 
y’’— psin= B = —}y8’ cosé. (21) 
The following discussion will be based on these two non-linear simultaneous differential 


equations for 8 and y.* 
In terms of 8 and y, again under the assumption that \ < 1, we have for direct and 


bending stresses and for displacements the following formulas 








pe See 
o> =", = Tea — vp bY? @) 
— Ny - aes Mn m [(cos— + Bsin dy + (sin — — B cos §)Q]) (23) 
m=) = ed — A) a hO ani ’ 
6M E i 

u 1 a. 

baa TY = LY - 
. = -|/ (4 cos —& + 56 sin :) dé. (26) 


4. Boundary conditions for problem of Omega joint. Since the solution of the problem 
in accordance with Fig. 1 is symmetrical with respect to the plane z = 0 (§ = +2/2) 
we need consider only the lower half of the toroidal shell (—1/2 < — < 2/2). At the 
outside circumference of the torus (—§ = 2/2), we have the symmetry conditions of no 


*Equations (21) and (21) may be considered as equations for the leading terms in a development 
in powers of \ of the solutions of the complete equations (17) and (18). 
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change of slope of meridian curve due to deformation and of vanishing horizontal stress 
resultant H. In terms of 8 and y these conditions read 


B(x/2) = 0, Y(r/2) = 0. (27) 
At the inner edge of the torus (§ = —2/2) we have conditions of elastic support 


through the two cylindrical shells to which the toroidal shell is joined. The complete 
formulation of these conditions is required when stresses and deformations are to be 
determined for small and moderate values of u. Our previous work [1] on the basis of 
linearised theory indicates that when yu > 1 the exact form of these boundary conditions 
is unimportant for the determination of significant stresses and deformations. 

Since explicit calculations in what follows are restricted to the case up > 1 the re- 
maining boundary conditions are taken in the form 


B(—7/2) = 0, y(—7/2) = 0. (28) 


5. Range of applicability of linearised theory for moderate values of parameter yu 
and method for determination of non-linear corrections. Inspection of the differential 
equations (20) and (21) indicates a difference in the character of the solutions depending 
on whether u is of order of magnitude unity or is large compared to unity. 

In the former case, that is when 

php = O()), (29) 


the only parameter which can become large is the load parameter uQ. This possibility, 
however, is found to be excluded by the physical conditions of the expansion joint 
problem. As long as (29) is satisfied we have the order of magnitude relations 


8B = O(uQ), y = O(Q). (30) 


From this the order of magnitude of the deflection w is 
w= O(buQ) (31) 


Since w = O(b), at most, it follows that uO = O(1) is the practically significant range. 

Introduction of (30) into (20) and (21) indicates further that the linearised theory 
is applicable as long as 

wQ K 1 (32) 

or, in view of (31), as long as the deflections are small compared with the radius b of 
the cross section of the torus, independent of the thickness h of the shell. 

Corrections to the linearised theory may be determined by developing the solutions 
of (20) and (21) in powers at the parameter 
(33) 


pQ = A. 
If we write 


B= wQ > AB, , y=2 > A'y, (34) 


n=0 n=0 


and substitute these expansions into equations (20) and (21) we obtain 


> (8 + siné y,)A" = cosé + A(sin g + >> y,A” cos :) > B.A", (35) 


n=0 n=0 n=0 


© @ 2 
> (Ww — p’siné B,) A” = -5 wa( 2.4") cos £. (36) 
n=0 
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Equating coefficients of like powers of A then gives the following pairs of equations for 


6, and y, : 
0 +siné ¥ = cosé, 


(37) 
y - u sin & By _ 0, 
n-1 
Bi’ +sint ¥, = Bi sine + 2) Ben-z-1 Cost, n>, 
k=0 
(38) 


. 1 n~1 
vi’ — vw sin=B, = 5 u° 2, BBs. cost, n> I, 
k=0 
where the first pair (37) are the equations for the linearised theory. 
For sufficiently small values of » each successive pair of equations (37) and (38) 


may be solved by developing §, and y, in powers of the parameter yu as 


B, = p> Bint, Va = > Vnmit”« (39) 
For moderate values of » another approach is to determine solutions of (37) and (38) 
as Fourier series in a manner similar to that used for Mathieu functions. 
As our work here is primarily intended to furnish an analysis valid when p > 1 we 
shall not carry out any explicit calculations for the case » = O(1). 
6. Approximate solutions for large values of 1. We now consider the system of differ- 
ential equations (20) and (21) with the boundary conditions (27) and (28) under the 


u> il. (40) 


Our approach is guided to some extent by what we know [1] concerning the solution of 
the linearised problem under the assumption y > 1. For the linearised problem we have 
a narrow zone around ¢ = 0 where bending is important. In the remainder of the toroidal 
shell bending is a secondary effect and membrane action dominates. We shall show that 
an approximate solution for the non-linear problem may be found which is qualitatively 


similar to the solution of the linearised problem. 
As long as we exclude the point § = 0 we may obtain formal solutions of (20) and 


assumption that 


(21) by expanding as follows 


B= > Bu", v= dD vu” (41) 


n=0 n=0 
Requiring the coefficients of all powers of » to vanish after substitution of (41) in (20) 
and (21) we have the following relations 

Bo(Bo — 2 tan &) = 0, 


(42) 

Yo = Q cote + Boy cotE + 08, , 

Ba. = Wil, esc § + 5 coté D Bibs; 
n> 1. (43) 


’ 


vn = —Bil, escé + 08, + cott D> Bina 
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From (42) we have, in view of the boundary conditions (27) and (28), 
8 = 0, vo = Qcoté, (44) 


which is the solution of linearised membrane theory. With (44) we have next from (43) 


B, = Wo’ esc é, y, QB, + WR: cot é (45) 
or 

B, = 22 esc* é cot é, ¥, = 22° ese’ — cot é. (46) 
In order to have the first two non-vanishing terms in the expansions both for y and 8 
we also determine 6, , with the result that 


0 eos 2/.. _ eos ¢ 
pa go oes E + (17 + 43 coef) a+ s+, 
usin & be sin’ é/ sin é 


Me £ oO 
} = 0 S088 | eee. 4 ma) 
sin € usin’ £ 


We may readily see that the complete expansions (47) satisfy the boundary conditions 
(27) and (28) term by term but their singularities at = 0 indicate that these expansions 
cannot be valid near é = 0. 

In order to obtain a solution which is valid near § = 0 we observe that for |~| <1 


(47) 


we may approximate (20) and (21) by 


BY + wey = w[Q(1 + BE) + YA], 


(48) 
vy’ — piB = —fus 
As for the linearised form of these equations we make the substitutions 
a= p's, (49) 
B= p'“Of(2), Y = pw” itg(z) (50) 
which transforms (50) into 
ff!’ +2ag=14+ ph Ofg + Azf, 
(51) 
g”’ — xf = — of’, 
where now primes indicate differentiation with respect to z. 
Te shall here restrict further attentions to problems for which 
We shall h trict furtk ttent to probl f hicl 
pw? = O(1). (52) 


Equation (52) evidently must hold in order that the solution of the non-linear system 
(51) should have the same qualitative behaviour as the linearised system. 
Since » > 1 we may now neglect the last term on the right of the first of equations 


(53). Setting further 
r = p09, (53) 
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we are left with the system 

f’ +29 =1+Tf9, g’’ — af = —1Tf’. (54) 
Since » > 1 we may replace the boundary conditions at = +: (x/2)y'”* by the following 
conditions 


f(+~) = 0, gio) = 0. (55) 


We assume the solutions of (54) may be written in the form 
f=Diior, g= 2D 9.(2)T". (56) 
n= n=0 
Substitution of the expansions (56) in (54) results in successive pairs of equations for 


the functions f, and g, . Each pair of equations may be combined into one complex 


differential equation for a complex function 7, defined by 
T(x) = f,(x) + t9,(2). (57) 


The resultant differential equations are 


Ti'(x) — 12T,(z) = R,(2), a: ©), 2. ++> (58) 
where 


R.(r) = 1. R(x) = > ag Lo. -- 5 if. ‘A n> i (59) 
A 1 


and where, in view of (55), 


T(t o) = 0. (60) 


Two linearly independant solutions of the homogeneous equation corresponding to 
(58) are given by the tabulated functions h,(7x) and h.(ix) which are defined in terms 


of Hankel functions [3] by 
kee re 
—* (Fe ‘ win(32), 


J ae. 
nae) = (2: ) A, (22 ). 


The method of variation of parameters furnishes the solution of (58), (59) and (60) 


(61) 


in the form 


in 9) 1/3 ar a & 
T.=- 4 (2) Ke | R,(t)h(it) dt + h.(iz) | R,(Oh, (ib) a], (62) 


where use has been made of the fact that the Wronskian of h,(z) and h.(z) is (4/mi)(3/2)'”. 

The function 7 , which is a modified Lommel-function, and its first derivative have 
been tabulated previously [1]. In order to obtain an idea of the magnitude of the non- 
linear effect we have calculated 7',(x) and its first derivative by numerical integration 
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using equation (62). Figures 3 and 4 contain graphs of the real and imaginary parts of 
the functions T, , 73, T; and T?. 

From the behaviour of the functions h, and h, it may be deduced that f,(z) is even 
in x when 7 is even and odd in x when n is odd while the reverse is true for g,(z). 

We note further that from the asymptotic behaviour of the functions h, and h, it 
can be concluded that the behaviour of 7, and 7’, for large values of z is as follows 


8 ! ” 
fo(x) ~ go(x) ~ ya (63) 
— 2 
ia~=, g(z)~ =. (64) 

1 x 


It remains now to interpolate between the functions 8 and y as defined by (50) and 
(56) for small values of ~ and the functions 8 and y as defined by (47) for values of & 
around § = +7/2. We shall carry out this interpolation for the first two terms in the 
series for 6 and for y. Write in equations (47), representing the solutions of (20) and (21) 
for sufficiently large values of &, 

y= p*sing (65) 


giving 


(66) 


bo 


, een 
y = u'*Qeosé Ll+-uep Qt ee I. 
y y 


We now compare this with the solution around = 0 which may be written, in view 


of (50), (53), and (56), in the form 


B= p?Olf(r) + wos (x) +--+], 
(67) 
Y= pw? O[go(x) + uw’? Qg,(x) + --°]. 
For sufficiently large values of = u'”£ we have from (63) and (64) 
pw eas, E + eS + | 
Lt i. 
(68) 


/ ] /3 2 
y~ n'a E + PO = + ++ 


x x 
In order that we may have a continuous transition between the solutions (66) and the 
solution (67) we now modify the solutions (67) to read 


B= pally) + wr’ Ofi(y) + --+], 
(69) 


II 


v= 2? OQ[goly) + uw’? OQg.(y) + ---]. 


Since z ~ y when £ is small we have agreement between (69) and (67) in this region. 





1953] A PROBLEM OF FINITE BINDING OF TOROIDAL SHELLS 331 





0.8}—- he 


0.6 F— 



































3.2 








332 R. A. CLARK AND E. REISSNER [Vol. X, No. 4 


~ 


For large y equations (69) are, in view of (63) and (64) 


(70) 


Equations (70) agree with equations (66) in all terms except one. The reason for the 
absence of one of the terms in (70) which occurs in (66) is due to the fact that we omitted 
the last term in the first of equations (51). While this is permissible as long as x is O(1) 
it is no longer permissible for | z | > 1 and so we cannot expect complete matching of 
the developments (68) around — = 0 and (67) around § = +7/2. Since as we shall see 
all significant numerical contributions of the solutions are to be found in the vicinity of 
§ = 0 this lack of complete matching may be admitted without effective consequences. 

7. Numerical results and determination of range of linearity for 1 > 1. In order to 
obtain an idea of the numerical effect of non-linearity for the expansion joint problem 
we must introduce the solutions (69) into the expressions (22) to (26) for stresses and 
displacements. 

Just as in the linear problem the significant stresses are the circumferential direct 
stress o¢p and the meridional bending stress o;g . These stresses have their largest values 


near the section ¢ 0 and in the vicinity of this section we obtain from (22), (24) and 
(69) 
E h oP ! 2/3 , j 7 ("1 
6 =F —3i72 FH Oi goly) te Qgily) + -°*], (71) 
121 —yv)]* b iia oe 
E h 2/3 “ \ 2/3q¢7 me 
Or; ; BQ oly) te QSi(y) + vee], (72) 
21 —yr)b Ts — 


where y ~ yu *E 

The corresponding results of the linearised theory are equations (71) and (72) with 
only the first terms inside the brackets [1]. 

Since gj and gj , and fj and f{ are of the same order of magnitude (see Figures 3 
and 4) we have in fact as the condition for negligible non-linear effects on stresses that 
the parameter »’’°2 must be sufficiently small. 

According to the linearised theory o¢p is an even function of — and has its maximum 
at & = O, while o;, is an odd function of £ and has its maximum for ~ = , 


where fé'(y,,) = 0 which relation is equivalent to 1 — y,,.go(y,) = 0. The value of y,, = 
u'“¢,, following from this relation was found in reference 1 to be y,, = +1.225. 
Inclusion of non-linear terms in (71) and (72) introduces both even and odd func- 
tions of ~ into agp and o;g . The maximum of o¢p then no longer occurs for y = 0. Let 
y, be the value of y for which o9p = o¢p, max. - Equation (71) indicates that 
sre 2/3 vr 
oY) te Og") ++: =O 

or, in view of the form of the differential equations for f and g, that 

¢ 2/3 - 2/3 me 

2y, — Kw Afoly) + wv’ Ofi(y.) + ---] = 0. (73) 


Numerical values of y, for a few values of u°’°2 may be found in Table 1, together with 
the corresponding values of g'(y,) = gé(y:) + w’’Ogi(y,). It is seen that consideration 
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of the first non-linear term increases the value of o9p, max aS compared with the value of 
the quantity according to the linearised theory. This increase is, however, only slightly 
more than two percent up to values of 0.3 for the quantity u’”°2. The location of the 
maximum of o»p is seen to be shifted toward the axis of the torus. 

Maximum values ¢¢g,max Of the bending stress o;g, occur according to (72) where 
fi (Ym) + WOOT (yn) + +++ = 0 or, in view of the differential equations (54) where 


Ym{Go(Ym) + we QGi(Yn) + +++] 
= 1+ pO folYm) + uw? OF (Ym) I[Go(Ym) + mw’? Qgi(y)] (74) 
L + we Q[fo(Ym)Go(Ym)] + --- 


Il 


Equation (74) has two solutions y,,,, and y,,2 . Values of y,, and y,,2 may also be 
found in Table 1. It is seen that with increasing values of u’**Q both of the locations of 
extremum values of o;, shift towards the axis of the torus. In order to see the magnitude 


of the non-linear effect on ¢¢8, max We have also tabulated values of f’(y,,.) = fi(Ym) + 
nu “Of {(y,,). Non-linearity is seen to increase the bending stress at the inside maximum 
location and to decrease it at the outside maximum location. When »?“2 = 0.3 the 


increase in stress is about six percent. When © is negative rather than positive this 


situation is reversed. 


TABLE 1: DATA FOR VALUES AND LOCATION OF MAximMuM STRESS 














p?!3Q ui g' (y:) Ym,1 f'(Ym,1) Um,2 Ff’ (Ym,2) 
0 0 0.939 —1.225 —0.753 1.225 0.753 
0.1 —0.064 0.942 —1.284 —0.767 1.155 0.744 
0.2 —0.130 0.949 —1.337 —0.784 1.080 0.741 
0.3 —0.200 0.969 —1.382 —0.804 1.000 0.743 





Let 6 be the relative axial displacement of the inner edges of the toroidal joint. We 


then have from (26) 


lL. T T = Dee fon) 
55 = wW5) —w “a —b | B cos—é + 9 8 sin £} dé. (49) 


v“=-9r/2 


With 6 from (69) and with y = u'” sin — equation (75) becomes 


. 1 
6=—-b| 4 y ay 


Snr dy | sacar 
pw OF(y) wat oH ‘Mf(y] 34 7 — (yp 


| 


= —2b | [(foly) + wOfi(y) + +--+) 


Q ee _ 
+ BY ai Go) + we Ofy) + +)" dy. (7) 


ene ee 
2{1 — (y/u 


Since f,(y) is even in y and f,(y) odd in y we see from (76) that 


§ = an | foly) dy + O(u* ‘ot |. (77) 
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The leading term in the omitted expression is due to the function f2(y) which has 
not yet been calculated. Since f; , which gave first order corrections for the stresses, 
does not contribute to the value of 5, we may conclude that non-linear corrections to 6 
are, for small values of u*’°Q, even less significant than are these corrections for the values 
of the stresses. We know [1] that 


| foly) dy ~ | foly) dy = —71 (78) 


and so 

6 = 2rQb[1 + O(u*”*0”)). (79) 
From (79) we may obtain the order of magnitude of the relative axial displacement up 
to which the linearised theory is applicable. The condition for linearity u’“Q2 « 1 may, 
from (79), be written in the form 


_ 2rb 2r (ab)*’° Be 
1S i (80) 


For practical purposes we may take instead of (80) 


9) I 
=~ a75- (80’) 


ou 


The inequalities (80) are remarkable for the facts that (1) displacements much larger 
than the wall thickness of the shell fall within the scope of the linearised theory, and 
(2) the ratio of permissible displacement 6 to the radius of the torus decreases as the 


parameter yu increases. 
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A BOUNDARY VALUE PROBLEM IN THE THEORY OF PLASTIC 
WAVE PROPAGATION* 


BY 
E. H. LEE 


Brown University 


Abstract. In the analysis of boundary value problems in the theory of plasticity, 
the general situation arises that different partial differential equations are to be satisfied 
depending on whether the material is in the plastic or elastic state. The criterion de- 
termining the state at any material point depends on the dependent variables and their 
derivatives with respect to time. Thus the regions of application of the different differ- 
ential equations must be determined from the boundary and initial conditions as the 
solution is developed. 

The theory of the propagation of plastic waves in one dimension is a case in which 
the solution, including the determination of the unknown plastic-elastic boundaries, 
can be treated. An example is presented in this paper which illustrates the many types 
of boundary determination conditions which must be used. The method is based on 
the numerical integration along the characteristics of the hyperbolic equations arising, 
one linear and one quasi-linear. The development is possible since forward integration 
along characteristics enables the unknown boundaries to be determined independently 
of the subsequent solution. This situation is contrasted with other problems in the 
theory of plasticity. The complexity of the procedure indicates the difficulty to be 
anticipated with analytical treatment of such problems, and with the numerical treat- 
ment of problems involving more extensive plastic flow. 

1. Introduction. The general mathematical theory of the propagation of plastic 
waves along rods was developed by Karman, Bohnenblust and Hyers [1].** The solution 
of boundary value problems in this field requires the treatment of a quasi-linear wave 
equation in the regions where plastic flow is taking place, but the form of the equation 
changes if the stress falls at any section which has previously undergone plastic flow. 
Such a region is termed an unloaded or hysterises region, and its existence and develop- 
ment depends on the solution of the boundary value problem prior to the instant con- 
sidered. This situation which requires the determination of the unknown boundaries 
between plastic and unloaded regions leads to a wide variety of boundary determination 
conditions. It is in this respect that the analysis differs from that of non-linear elastic 
waves in which the same differential equation would arise for both loading and un- 
loading, so that the unknown boundary problem would not appear. The solution de- 
tailed below is presented since, in view of the wide variety of boundary determination 
conditions needed, it illustrates in a general way the types of problems to be faced in 
this field, emphasising the difficulty to be anticipated in the analytical attack on such 
problems. 

The method of solution is based on the theory of characteristics of the hyperbolic 
equations arising, and is treated numerically. It was found that in order to maintain 
sufficient accuracy to determine the unknown boundaries between plastic and unloaded 
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regions it was necessary to compute the geometry of the slip line field using six sig- 
nificant figures. Graphical methods were found to be not sufficiently precise to de- 
termine consistent boundaries. Thus this is an example of the situation which may 
arise in the numerical treatment of boundary value problems, where more accuracy is 
needed to carry through the method than is warranted by the initial physical data: in 
this case the stress-strain relation of the material under consideration. It presumably 
means that in an experimental check minor variations in the material would involve 
appreciable changes in these boundaries. However in this case only minor changes in 
the final strain distribution would be expected to result. 

Because of the accuracy required to complete the solution, it would be very lengthy 
to treat a problem involving more extensive plastic flow. Moreover, because of the 
variety of boundary conditions encountered, such problems do not lend themselves to 
routine computing. 

2. Basic Theory. The equations governing the propagation of plastic waves were 
developed by Taylor [2] and Karman [3] independently. Karman’s theory has recently 
been published with some experimental results (Karman and Duwez [4]). The theory 
was later further developed principally by Bohnenblust in several O.S.R.D. reports. 
Here we shall consider the final form of the theory which applies for finite strain. 

We are concerned with a material having an invariant relation between stress and 
strain under continued loading, independent of the speed of the test. Such a relation 
is shown by OABD in Fig. 1. In the elastic region OA the relation is linear with gradient 


Of B . ? 








FIG. | 


of magnitude Young’s Modulus E. The are ABD represents plastic flow. The stress 
is defined as the nominal stress, force per unit initial cross-sectional area, and ¢ the 
nominal strain, change in length per unit initial length. These nominal values are used, 
since in combination with the Lagrange type space coordinate z, large deformations can 
be analysed without introducing additional complexities into the equations. For con- 
venience with the particular problem to be considered compressive stress and strain will 
be taken as positive. If when the point B has been reached unloading takes place, elastic 
strain increments only occur, and the stress-strain point moves on the straight line BC 
in Fig. 1 which is parallel to the original elastic line OA. If the stress later exceeds the 
value at B, the stress-strain point continues along the plastic flow curve BD. 
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Longitudinal motion in the z direction only is considered, and this puts a limitation 
on applications if the associated lateral motion becomes important. The equation of 


motion is 


Oa ou 
cy oa 2) (1) 
Ox ot 
where p is the initial density, and u the particle displacement from its initial position 
xz in the unstrained rod. If straining is taking place on the curve OABD, o = a,(e), 
¢ = —du/dx; d0/dx can be replaced by —o/(e)d’u/dx*, and the quasi-linear wave 
equation is obtained: 
Ou Ou a i do, ( os) 
—-_ =o —-; c= = f(e-) = fi ——}]. (2 
ot Ox” p de I P Ox 
The characteristics of this equation (see [1]) are given by dz/dt = Fc, and along 


them we have the integral relation 


dx 


= Fe, v + ¢ = constant, (3) 
( 


where v is written for the particle velocity du/dt, and g(e) = fic de. 
When the material has been unloaded, the corresponding equations have been given 
by Bohnenblust [5]. If B in Fig. 1 is given-by (€max , max), the unloading line BC is 


oo 7, (Emax) is E (nes a €). (4) 
The equation of motion becomes 
ou » Ou d|1 2 s 
at == ax’ a dx 5 70(Emax) — Co€max |; (5) 
where c, = E/p. The characteristics, and characteristic relations are 
dx : = 
= Fey, pc)’ = o = constant. (6) 
tf 


At any specified time, the second term on the right hand side of (5) is a function of 
x which is not prescribed initially, but which is determined from the solution at earlier 
times. 

We have summarised the two types of equations to be considered. They can be 
applied directly to continuous stress fields as long as the stress strain relation OABD is 
concave downward. When this is not the case it has been shown by White and Griffis 
[6] and Lee [7] that shock waves are predicted, but we shall not be concerned with this 
situation in the present paper. However, in the elastic and unloading regions waves of 
stress and strain discontinuity can be propagated along the characteristics, and the 
Hugoniot relations of continuity and momentum change take on forms very similar to 
the characteristic relations. We shall be concerned with such finite linear waves. 

We shall consider the solution in terms of the characteristic net in the (z, ¢) plane 
and for convenience we replace the time coordinate by r = ec ot/l, and xz by —& = 2/l, 
where | is a typical linear dimension. The linear characteristics are then given by 
dt/dr = #1. The characteristic relations in the unloaded region U are then 


dé 


= ¥1, pC’ - o = const. (6a) 
dr 
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and in the plastic region P 


a = Fc’, pCov - peop = const., (3a) 
where c’ = ¢/¢y. 

3. The impact problem considered. The problem we shall consider in detail, which 
involves an illustration of the many types of free boundary determination conditions 
which can arise in the analysis of plastic wave boundary value problems, is the motion 
resulting from the normal impact of a cylinder of length 1, moving parallel to its axis 
with velocity v) , against a rigid target at rest. For the convenience of having zero initial 
conditions we add a constant velocity v) to the whole system. The initial and boundary 


conditions are then 


= 0, og = U, T= Q, 0< g <ul, 
v=, §€=0, 7>0, (7) 
c= 0, é = . T > QO. 


9 


The initial development of the corresponding characteristics field is shown in Fig. 2. 
The discontinuity in velocity at 0 combined with the zero initial conditions determines 
a simple wave or lost solution with a fan of straight characteristics through 0 as dis- 
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FIG.2 
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cussed by K4érman and Duwez [4]. Along each straight characteristic the velocity and 
stress remain. constant. OA is the elastic wave front across which the stress rises in- 
stantaneously from zero to the yield stress Y, corresponding to A in Fig. 1. OD is the 
characteristic corresponding to the velocity v) . Throughout the triangle ODE the ma- 
terial is moving with the constant velocity v, and is subjected to constant stress. 

At the free end, = 1, the elastic loading wave is reflected at A as an unloading 
wave of stress discontinuity ABC, which is propagated with the elastic wave speed Co 
in view of the unloading relation, and so has gradient —1 in Fig. 2. At A the reflected 
wave is equal in magnitude and opposite in sign to the original elastic wave front and 
produces zero stress there, but as it advances into the plastic region its magnitude is 
reduced as discussed below. Let the suffixes b and a denote conditions before and after 
the unloading wave ABC has traversed an element of the cylinder. Below AC the 
material is in the plastic condition P, while above AC unloading has taken place and 
the material is in the condition U. At any point along AC o, , «, and », are known from 
the simple wave solution. The momentum relation for the unloading shock wave corre- 
sponding to the characteristic relation (6) is 


eo %& = — plo(V, _ Up). (8) 


At A pctov, = o, = Y for the elastic part of the initial wave solution. Since ¢, = 0 at A, 
(8) determines pcyv, = 2Y there. Considering the characteristic adjacent to AC, the 
constant in (6a) is given in terms of the known conditions at A, and we have 





PCV. — Og = 2Y. (9) 
Eliminating v, between (8) and (9) gives 
i<a« ace (10) 


Application of (3) for the initial wave solution emanating from 0 gives 
nee re 
PCV, = Co | cde = | — de. * (11) 
Jo Jo € 

since do/de = pc’. But c = Cp in the elastic region OA in Fig. 1, and c < ¢» when plastic 
flow occurs, and so it follows from (11) that the quantity [pcov, — o»] in (10) is zero at 
A in Fig. 2 and increases from A towards C. Thus (10) shows that the unloading shock 
wave is progressively absorbed as it penetrates the initial plastic wave front. If the im- 
pact velocity is below a certain critical value determined by (10), the unloading wave 
of stress discontinuity will continue through to the impact fact at £, but for higher 
impact velocities the unloading wave will be absorbed by the original plastic wave 
front. Thus, the stress discontinuity across the unloading wave may be reduced to zero 
at the point C, say. 

Thus two types of problems arise in analysing the subsequent motion. The plastic 
region may spread from C above AE, or the material may be unloaded throughout its 
length above AC. In the latter case it does not follow that plastic flow does not occur 
again. If at any section the solution of the equations for unloaded motion determines a 
stress equal to the previous maximum stress to which the section has been subjected, 
plastic flow can occur again there, and a new region in which the plastic wave equations 
must be solved is initiated. Such a problem is illustrated in this paper, and the dis- 
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cussion of the determination of the boundaries of such a region is given in the next 
section. 

To illustrate the results in detail, the solution is presented for impact of a cylinder 
of aluminum having the stress-strain curve shown in Fig. 3. The corresponding ¢, ¢ 
relationship [see Eq. (3)] is shown in Fig. 4. For this material the unloading wave of 
stress discontinuity is propagated through to the impact fact for impact velocities 
vo < 363 ins/sec. For values greater than this the unloading wave is absorbed and the 
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plastic region spreads out above it. The case considered in this paper, the characteristic 
field for which is illustrated in Fig. 5, corresponds to an impact velocity of 356 ins/sec 
so that the unloading wave just gets through without being absorbed, but for which 
plastic flow occurs again subsequently. The unloading wave is reduced from the yield 
point stress magnitude, 5 X 10° p.s.i., at A in Fig. 5 to 0.131 X 10° p.s.i. along BC 
where it remains constant. 

4. Determination of the plastic-unloading boundary. For a material with a stress- 
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strain curve which is concave to the strain axis as considered here, shock waves of stress 
discontinuity are not initiated within the field of solution but must emanate from dis- 
continuities in the boundary or initial conditions. Thus, if the unloading shock wave 
ABC in Fig. 2 has been absorbed at C, the subsequent solution of the impact problem 
specified in (7) does not involve waves of stress discontinuity. However, discontinuities 
in stress and velocity derivatives may occur, and in fact are in particular associated with 
the continuation of the plastic-unloading boundary from C. Such discontinuities must 
therefore be considered in investigating the conditions which determine the boundary. 
A similar situation arises when a new plastic region is formed if the unloading wave 
ACE in Fig. 2 traverses the entire cylinder. The conditions determining the plastic 
unloading boundary have been developed by Karman, Bohnenblust and Hyers [1], but 
since this report may not be readily accessible, a brief statement of the theory is given 
in the Appendix with a somewhat simplified derivation. 

In the previous section we considered the characteristic field shown in Fig. 5 up 
to the unloading shock wave AC. The solution beyond this can formally be continued 
by assuming the maintenance of the unloaded condition and making use of the char- 
acteristic relations (6). Before this extension can be considered satisfactory, it is neces- 
sary to check that the stress at any section does not exceed the previous maximum value 
reached just ahead of the unloading shock wave AC. This is, in fact, found to occur 
at D. The characteristic OB determines the wave front of maximum stress, behind 
which in OBC the stress is constant and equal to the maximum impaci stress. In the 
unloaded region the influence of the wave front and the discontinuity of stress gradient 
across OB is transmitted along the characteristic BD, leading to the initiation of the 
new plastic region at D. A similar condition for which according to the unloaded solution 
the stress reaches the new vield stress at each section determines the plastic-unloading 
boundary DE. This is a boundary of type (i) as defined in the Appendix, and since 
dx/dt > cy (d§/dr > 1) no other conditions are required, and the two intersecting char- 
acteristics in U and the known value of o,,,. = o, on AC determine its position. If these 
conditions are formally extended beyond £, it is found that dr/dt < co , and a new 
condition discussed in the Appendix must be used. It is, of course, evident that the old 
cannot be applied in this region since the characteristic dé dr = +1 intersects 
» can no longer be used to determined the boundary. However, 
the formal continuation of the boundary DE beyond E is useful to obtain the termination 
point E accurately. Since the boundary at F is tangent to a characteristic d§/dr = +1, 
+ g) is zero there, where s is the are length along the boundary, so that the 


conditior 
the plastic region and 


0/08( plot 
position of E can be determined from a plot of (pcov + o) extending on both sides of E£. 
From the known conditions along the boundary DE the solution in the plastic region 
is obtained by the method of characteristics for a Cauchy problem using finite differences 
and the characteristic relations (3). 

A power series expansion of the variables about E indicates that there is a discon- 
tinuity in the gradient of the boundary there. In fact dx/dt becomes less than c, and the 
boundary EF is determined by a characteristic in each of P and U, and the known 
value of the yield stress at each section om; . The determination of EF permits the 
plastic field to be extended. At F the boundary becomes vertical, so that a continuation 
would involve a type (ii) boundary as discussed in the Appendix and modified conditions 
must be satisfied. It is found that dz/dt becomes <c, and the condition dr/dit = 
—1/p (do/dv) must be used with a characteristic in each of P and U to determine FG. 
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The upper part of the boundary DK of the plastic region P emanating from D is 
determined by two characteristics in P and one in U. dx/dt > c, so that with this type 
(ii) boundary no further conditions are to be satisfied. At K this boundary intersects 
the remnant of the unloading wave reflected from the impact surface at C. The P-U 
boundary then takes the form of an unloading wave of stress discontinuity with dz/dt = 
¢) along KJ. The known solution in P, combined with the momentum relation and the 
characteristic value on the U characteristic KJ determines the magnitude of the dis- 
continuity across the unloading shock, and this is found to be completely absorbed at 
J. Thereafter, along JG, dx/dt < c) , and two characteristics in P and one in U deter- 
mine the segment of the boundary which closes the plastic region. 

A check of the characteristic constants determined by the solution integrated to the 
point G, shows that the new yield stress at each section is not exceeded in the subsequent 
motion, so that the U region covers the entire cylinder for all later terms. The cylinder 
will rebound in elastic vibration, with a permanent plastic strain distribution generated 
by the two plastic regions P. 

5. Discussion of the solution. The permanent plastic strain distribution corresponding 
to the characteristic field of Fig. 5 is shown in Fig. 6. It is obtained from the maximum 
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strain reached at each section when plastic flow finally ceases there, by subtracting the 
elastic component of strain o,,./H. ABCGEF is the plastic strain produced by the 
impact wave spreading from 0 in Fig. 5 before it is absorbed by the unloading elastic 
wave reflected from the free end. The subsequent flow due to the second plastic region 
produces the additional permanent strain CGED. This is seen to be small compared 
with that due to the primary plastic wave. The major part of the strain is seen to be 
concentrated at the impact end, which for higher velocities and large strains leads to 
the mushrooming effect of impact. 

A study of this solution indicates the extensive detailed analysis required to de- 
termine the boundaries of the secondary plastic region even though it is quite limited 
in extent. It indicates the lengthy calculation to be expected for a case of higher velocity 
impact where detailed boundary considerations will be necessary for a much more 
extensive plastic region. The changing type of boundary, with the different conditions 
discussed in the Appendix, prevent a routine computational program and so increases 
the computing time. High numerical accuracy must be maintained in order to determine 


consistent boundary points, particularly on the segments where the condition dx/dt 
—1/p (do/dv) is used, since it involves, in its finite difference approximation, operating 
with the small difference between almost equa! numbers. In the present case it was 
found necessary to work with six significant figures, and to compute all characteristic 
intersections in the z-¢ plane, rather than making use of a graphical plot. 

This problem illustrates the essential difficulty which arises in boundary value 
problems of the theory of plasticity. The different relations governing elastic deformation 
and plastic flow lead to different differential equations for material in these two states, 
and the domains through which each of the two states pertain must be determined from 
the prescribed boundary values on the outer boundary. This poses a much more difficult 
problem than the usual boundary value problem for a prescribed differential equation. 
The present problem of plastic waves in one dimension is a case that can be handled for 
general boundary conditions by the methods discussed above. The solution is made 
possible by the hyperbolic nature of the wave equations governing each type of be- 
haviour, which permits forward integration along real characteristics, and the de- 
termination of the domains in which each equation applies at each stage as the solution 
proceeds with increasing time ¢. In other problems in the theory of plasticity such a 
direct procedure is not possible. For example, in plastic-rigid theory of quasi-static flow 
in plane strain, which has received considerable attention in recent years, different 
conditions are to be satisfied in the plastic and rigid regions. In the plastic region two 
hyperbolic pairs of first order partial differential equations arise which have the same 
pairs of characteristics, the slip lines. In general the prescribed boundary conditions are 
such that insufficient boundary information is available for forward integration of both 
pairs of equations, so that the solution at some point near the boundary may depend 
on the plastic-rigid boundaries to be determined within the domain. Thus, in general 
only very special problems can be treated, in which for example static determinacy 
permits separate integration of the stress pair of equations from the prescribed boundary 
conditions. From the point of view of considering the general types of boundary value 
problems in the theory of plasticity, it is useful to have one group of problems, waves 
in one dimension, for which general boundary value problems can be treated directly 
as illustrated in the present paper. 

Acknowledgement. The solution discussed above was taken from an official report 
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Appendix. Determination of the P-U boundary. 

As discussed in Section 4, we are concerned with a boundary across which derivatives 
of stress and velocity are discontinuous. For brevity we will use subscripts 2 and ¢ to 
denote partial differentiation with respect to these variables, and superscripts p and u 
to denote values in the plastic and unloading regions. 

Continuity implies that the differentials must be the same on the two sides of the 


boundary, therefore 


do =crdr+oardt = oa,dr+a; di, 
(a.1) 
dv = vdxrt+ vidt = vidxt+ vy; dl. 
Making use of the equations of motion and continuity ¢o, = —pv, , 0, = —pcv., 
all derivatives but o, can be eliminated from (a.1) giving: 
\o 
1 (/dx\~ 
pT a Nat 
adi = = (a.2 


— = ee ee 
oc 4] (42) 
c \dt 
0 in order to maintain plastic flow. Two cases arise: (i) an element on 


the boundary is changing from U to P, (ii) from P to U. They are considered below 


in which oa? : 


side by side. 


i UP, (ii) P— U, 
o, > @: a = ©. 
From (a.2) a > 0, From (a.2) a < 0, 
dx dz 
ae or <ec Cc — st ee 
dt 3 = ee 
otherwise o? = ao; = 0. otherwise o? = o; = 0. 


When oa? = o; = 0, it follows from (a.1) that either ¢, = 0 or 
ldo _ _ dx 
pd at’ 


These conditions serve to determine the boundary as the solution is evaluated by 
characteristic integration. 
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MINIMUM WEIGHT DESIGN AND THE THEORY OF PLASTIC COLLAPSE* 
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Summary. This paper examines the problem of assigning economical sections to 
the members of a structure whose geometrical form is given. The criterion of failure is 
taken to be that of the plastic theory of collapse, and the criterion of minimum weight 
is employed to determine the best design. A geometrical analogue of the equations in- 
volved is used to clarify their significance, and such proofs as there are in the text, are 
cast into geometrical terms. A method of solution is suggested at the end of the paper, 
but the primary concerns of the paper are the general features of the problem. 

1. Introduction. The problem which this paper considers may be stated as follows: 
given the geometrical form of a plane frame, together with a set of static loads which 
bear upon the frame, what cross sectional dimensions must the members be given in 
order that the loads may be just sustained and the weight of the structure as a whole 
be as small as possible? The members are assumed to have constant cross sections 
throughout their length. 

In order to decide whether or not a frame will just sustain the loads, a criterion of 
strength must be assumed, and the criterion employed here is that of the theory of 
plastic collapse [1].** This theory is applicable to framed structures of ductile material 
which derive their strength from a bending action, e. g. rectangular building frames, 
pitched roof portals etc. It assumes that the bending moment in a member tends to a 
maximum limiting value as the curvature becomes indefinitely large. This limiting 
value of the bending moment is called the fully plastic moment, and if this moment is 
attained at a particular cross section a plastic hinge is said to have formed there. A 
plastic hinge allows a finite change of slope to occur at the cross section where it forms, 
and the fully plastic moment resists any further rotation of the hinge. 

If the collapse of a structure is now defined as the condition when deformations 
increase indefinitely under constant loads, then, assuming that instability effects can 
be ignored, it can be shown that collapse occurs when a sufficient number of plastic 
hinges have formed to transform structure, or any part of it, into a mechanism. A load 
factor \ is usually employed to determine the correct collapse mechanism. Suppose, for 
a frame, whose members 7 have sections such that their fully plastic 
moments are §; , is subjected to the loads W, acting at the points k. (The notation used 
omits a summation symbol; all small suffixes should be read as repeating ones, capital 
suffixes refer to particular variables or equations). Then, if an arbitrary mechanism 
Y is assumed, an equation of virtual work may be written down 


AyWid&y = B96; Y 


where 6,» are the displacements of the loads and 6;y are the hinge displacements of the 
members, relevant to the mechanism Y. This equation determines the value of Ay and 


instance, that 
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it can be proved that Ay > Ac where X¢ is the load factor for the correct collapse mechan- 
‘ism C of the frame. [2]. It can be seen that if the frame is to withstand the loads W, 
the 8; must be such that Ac > 1 and thus Ay > 1. This briefly is the concept of strength 
used in this paper, so the assumptions above are inherent in all that follows. 

J. Heyman has already suggested a method of solving specific problems of this type 
[3] but tl 
it with t 
Economics [4]. 

The text is divided into three sections. The first states the general problem and the 


1e purpose of this paper is to clarify the nature of the problem, and to identify 
he problem of “activity analysis’ or “linear programming” which occurs in 


second introduces a geometrical analogue of the equations involved which helps to 
clarify their significance. The last section suggests a method of solution, but this is 
included for the sake of completeness, it is not recommended as a practical design pro- 
cedure. 

I. The problem. Suppose it is required to design for minimum weight, a frame 
whose members 7 have lengths /; when it is subjected to the loads W’, acting at the points 
k. If sections are assigned so that the fully plastic moments of the members are £; it 
will be assumed that the weight, G, of the frame as a whole may be represented by the 
expression 

G = 1,8 
iB; 

The assumption here, of course, is that the cross sectional area A of a member is 
proportional to its fully plastic moment. In fact A « 6" where 2 > n > 1 for most 
sections, but the assumption ‘that n = 1 -does not seriously restrict the theory which 
follows and its implications will be pointed out at the end of the paper. 

Now suppose that with these values of 6; a collapse mechanism Y is assumed for the 
frame. The equation of virtual work for this mechanism is 


and if the frame is not to fail in the mode Y, 


This condition may be rewritten, 


ay;b; > i, where ay; = 


and, since the equation is one of virtual work, all terms are positive, and so all the a’s 
are positive. 

This equation is in effect a constraining condition on the possible values of the 8, 
and similarly every conceivable mechanism of the frame will yield one such constraining 
condition. Hence, if all conceivable mechanisms form a set of conditions j the constraints 
of the problem can be expressed by the sets of inequalities. 
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when the variables 6; are subject to the restriction 
a; :B; > l 


all terms in the expressions being positive or zero. 

Considered as a problem in linear algebra, the problem is difficult for the reason 
that the significance of the constraining inequalities is difficult to grasp, but if the problem 
is translated into geometrical language their meaning becomes clear. 

2. The analogue. The basis of this geometrical analogue is to give the variables 8, 
the significance of coordinate variables in an n-dimensional space. But before proceed- 
ing to the general case it would be as well, perhaps, to examine a simple two variable 
problem from this point of view. 

Example. Let it be required to minimise the weight of the continuous beam shown 
in Fig. 1(a) when it has to bear the loads shown. The loads, and the lengths of the spans, 
in the problem have been divided by convenient factors in order to obtain simple numbers. 

Let the left hand beam have a fully plastic moment of 8, and the right hand beam 


one ot £2. 
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In this example there are four possible hinge positions, and there is one redundancy. 
Hence there are three basic mechanisms [5]. These are sketched in fig. 1, and it is easily 
seen that six different mechanisms can be obtained by superposing these three in various 
ways. There are thus six conceivable mechanisms in this example and the load factor 
expressions for these mechanisms are as follows: 
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Now if 8, and 8, are regarded as coordinate variables, the various load factor ex- 
pressions: 
a8; + a8. = A, 


will form sets of lines whose perpendicular distances from the origin will be proportional 
to \, and if the \, are set equal to unity, Fig. 1(e) is the result. Now since all load factors 
must be greater or equal to unity, the values of 8, and 8, which will allow the loads to 
be carried are given by points which lie above and to the right of all these load factor 
lines, in the region marked ‘‘permissible’”’. This permissible region is outlined by the 
load factor lines which correspond to the most critical mechanisms for the beam, and 
any load factor line which does not touch this region, e. g. 4% and Aé corresponds to an 
impossible mechanism, a mechanism which, regardless of what values 6, and , are 
given, will never occur. The reason for this is that modes such as \% always violate the 
yield conditions, or, in mathematical terms, the inequality \% > 1 is dominated by other 
inequalities such as Ac > 1. 
The expression for the weight of any particular design, 


G = 26, + 282 
is again a line on the diagram whose perpendicular distance from the origin is proportional 


G and so the comparative weights of structures given by various points on the diagram, 
can be judged by their perpendicular distance from the line through the origin 


28, + 26, = 0. 


It is clear then that a minimal beam is given by the point, or points, of tangencv 
when the line 
G = 28, + 28, 


is tangent to the permissible region. In this case the single point M is this point of tan- 
gency, and the values of 8, and B, at this point are 


B, = 13, B. = 3. 


It is possible in some cases to have a line of tangency instead of a single point. This 
happens in the example of fig. 2 and, in this instance a range of values of 8, and 8, will 
give minimal frames. 

It can be seen from the above figures that the permissible region is purely convex 
toward the tangent weight line, there are no reentrant angles to the region, and hence 
there is only one minimal position for the weight line. 

3. General theory. In general if the 8; are regarded as coordinate variables in an 
n-dimensional space, the equations 


a; 8B; = 1 


will form a set of flats which will divide the positive section into a number of regions. 
If the external side of a flat is defined as that side which does not contain the origin, 
then a region can be defined, called the “permissible subspace,” which is external to all 
flats. This subspace contains all the points which assign values to the 8; such that 


a; Bi > 1, 
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i. e. the points within or on the boundary of the permissible subspace represent struc- 
tures which will bear the given loads. 

The following theorem lays severe restrictions on the possible forms of the permissible 
subspace. 

Theorem. If 8, and 8, are any two points within or on the boundary of the permissible 
subspace, then the straight line joining these two points lies wholly within the per- 
missible subspace. 

Proof. 8, and 8, are external to all flats. If in proceeding from 8, to 8, in a straight 
line a flat X is crossed, a region is entered which is not external to flat X¥. Thus X must 
be recrossed at least once more before the line proceeds to £, , for 8, is external to flat 
X. A straight line, however, cannot cross a flat more than once, and thus the straight 
line 8, 8, cannot cross any flat, and so the theorem is true. 

The statement embodied in this theorem is the mathematical definition of convexity, 
and so, since 8; > 0 and a;; > 0 it can be shown that the permissible subspace is a convex 
set contained in the positive sector. Further, by employing the fact that a > 1; > 0 
it can be proved that there is always one, and only one, position of the weight flat in 
which it is tangent to the permissible subspace. Hence there is always one, and only 
one, minimal solution to a given problem. The solution may of course allow a range of 
values for the 8; , in which case all the designs of the range have the same weight, and 


all are minimal solutions. 
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4. A method of solution. It would seem that any method of solution must involve 
either a total, or a partial, guided exploration of the field of the positive sector in the 
above analogue. The former mode of attack is the one employed by J. Heyman in his 
inequalities method, and the economists also favour this line of attack. Any such method, 
however, suffers from the serious disadvantage that all the constraining inequalities 
a,; B; > 1, i.e., all the conceivable collapse mechanisms of a frame, have to be accounted 
for in algebraic terms. E.g. In the example of Fig. 1(a) all the conceivable mechanisms 
could be accounted for by the formula 


6*(3B, —— 3) + 64(38. aes 1) + 64(6,; a? Bo) > QO, 


where the 6*’s are Kronecker deltas which can have values of one or zero. If all the 
possible combinations 6 = 0 or 1 are taken, all the constraining inequalities will be 
accounted for. The task is thus not impossible, but prohibitively laborious in any but 
the most simple cases. 

The method to be described here involves a partial exploration of the field, which 
discovers the constraining inequalities as it requires them. 

Suppose Fig. 3 represents in diagrammatic form a section of the permissible region, 
and suppose point P represents an assumed form of the frame 6% . If the weight of this 
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form is Gp and the collapse mechanism is Y, P is a point on the intersection of the two 
flats 
G = 18%, Ay = ay,B% = 1 
Now if 8% is given a small displacement 6 8% in the flat Gp in such a sense that | Ay | 
is increased, P moves along Gp in the figure to the point Q. A recalculation of the collapse 


of 











(9) 
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mechanism and load factor will then move Q down to R. If 6 8* is too large, however, 
the path P Q’ S may result. But, in this case, a further operation S T R’, will still result 
in a lowering of the weight flat from Gp to Gz . If the operation P Q R or the whole opera- 
tion PQ S, S T R’ is referred to as a step, the method consists in making a series of such 
random steps until the minimum is found. 

Example. It is required to design the frame of Fig. 4(a) when it is subjected to the 
loads shown. Let the ringed numbers designate the members of the portal. 

Assume as a starting point that the sections of the members are in the proportions 


6, = 1, 8, = 1, B; = 1 


The collapse mechanism can be determined as that of Fig. 4(b), [5], and for unity load 
factor the sections are 


This is the point P referred to above, and it lies on the intersection of the two planes 
Gp = 28,+ 82+ 83, 
Ay = 38, + $62 + $6; = 1 


Now a displacement 6 8* is required which will keep Gp constant and will increase | \y |. 
There are various ways in which this can be done and what these various procedures 
imply can be seen in Fig. 5. Any displacement takes P to P’, PP’ lying in the Gp plane. 





Fic. 5. 
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The particular displacement chosen will specify the angle @ that PP’ makes with the 
line of intersection of the planes UV. The method of giving displacements could be 
systematised by specifying that PP’ had to be at right angles to UV as PP? is drawn, 
and a succession of such defined steps will usually be fewer in number than a series of 
random ones. Also, by further restricting the nature of a step the process can be guaran- 
teed to be a converging one. But, the extra calculation required to ensure a particular 
type of step invariably consumes much more time than the extra steps necessary in a 
random search. A series of random steps is thus recommended, and the example will 
be solved by such a method. 


If we write the weight and load factor expressions as follows: 


G = 18. , Ay = ay;B; 
where Y is the relevant mechanism at any stage, then the only rules required in making 


a random step are: 
values 8 with a large quotient ay/] are to be increased, those with a small one 
are to be decreased. 
it) The increments added or subtracted to or from the values 8 are to be in the 
ratio of the lengths / so that the weight is kept constant. 
It is also advisable not to alter more than two variables in any step. 
Proceeding by this method then, it is patent that 8, has to be increased and 8, de- 


creased. Hence, 


Step 1. Increase 8, by 3 decrease 8, by 3 then 
5 t 2 
5 — » 9) 8 > > 
o oO ) 
A check shows that the mechanism is now that of Fig. 4 c) which gives the load 
factor expression as 
l - 
A J 9 _ Dz l 
} 8 8 
Accordingly, the value > remain as they are. 
Step 2. The above load factor expression indicates that 6. has to be increased and 
8, decreased. Therefore, increase 6, by 3 decrease 8, by } then 
5 5 l 
3, = 35; 3. = 5; B3 =>5 
3 3 3 
These values give the mechanism as that of Fig. 4(d) which has the load factor ex- 
pression 
2 Se a i ] 16 
A = ae nm = - 6, = py 
o 5 o 15 
and hence ior unity load factor 
25 25 5 
3 = & —_ B3 = 
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It is now seen that the coefficients of the load factor expression are in the same ratio 


to one another as those in the weight expression: 


G = 28; + B» + B3 , A = 


grid 


l = 
B.+eh+=Bs. 
5 5 


No further steps are thus possible and the operations have reached the minimal 
solution. This has occurred in this particular problem because there existed a load factor 
plane, 

2 ] I 
A=78,.+7-h+-8 = 1 
5 5 5 

which was parallel to the weight plane, and was not an impossible one. Clearly in the 
tangential position, the weight flat will always adjoin those possible load factor flats 
which are least inclined to the weight flat, and if there exists a possible load factor flat 
which is parallel to the weight flat, the minimal solution is the range of values of 8; 
which lie on this flat within the permissible subspace. The interesting conclusion may 
be drawn from these facts that for a frame with a particular geometry there are certain 
mechanisms which the minimal solution will favour, e.g. In the case of a portal with 
the relative proportions of that above, the mechanisms of Figs. 4(d) and (e) give load 
factor flats which are parallel to the weight flat. So, if the loads will “allow” either of 
these two mechanisms, that mechanism will be the minimal collapse mode. If both of 
these are impossible, the minimal solution will probably be the line of intersection of 
the planes relevant to mechanisms of the type shown in Fig. 4(c) and (6), i.e. the section 
of the beam of the portal would equal the section of the right hand stanchion. The 
plastic hinge required at the upper right hand joint may therefore form either in the 
beam or in the stanchion, that is, the frame has two alternative modes of failure. In 
general if there are n independent variables 8, , a minimal solution will be a corner of 
the permissible subspace where n flats meet. A minimal solution, not necessarily the 
only solution, will therefore have n alternative mechanisms of collapse. This implies that 
the members meeting at a joint must be proportioned so that the joint can fail in either 
of two ways, as above, and/or that some members have to be proportioned so that their 
failure or non-failure makes no difference to the value of the load factor. 

Conclusion. This paper, it is hoped, has clarified to some extent the nature of the 
problem of minimum weight design, when the theory of such design is based on the 
assumptions of the theory of plastic collapse, and on the assumption that the weight 
of a member is proportional to the product of its length and fully plastic moment. This 
last assumption does not seriously limit the validity of the above theory, for the true 
expression for the weight of a member merely converts the load factor flats above, into 
slightly curved surfaces; the general features of the problem are not altered at all. 

A more serious limitation on the above theory is the assumption, taken over from 
the theory of plastic collapse, that instability effects can be ignored. For it sometimes 
happens that the loading conditions of a frame will force the theory to choose long 
slender columns which will obviously buckle. This occurs when a rectangular frame has 
little or no sidesway loads, in ordinary cases quite reasonable sections are usually de- 
manded. This danger can be offset, however, by stipulating certain minimum values 
for the 8; of the stanchions and setting the theory to work with these additional re- 
straining conditions. This is easily done for instance in the method of random steps 
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above. Another advantage of the method is that the calculator can stop his calculations 
at any stage and still have a design which is adequately strong and lighter than his original 
guess. The drawback of the method is, of course, the laborious and uncertain nature of 
the calculations. However, the calculations involved in any other method to date have 


t 
r 


o be programmed for calculating machines, and there is no reason why the method of 
andom steps could not also be similarly programmed. 
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AN EXTENSION OF THE METHOD OF TREFFTZ FOR FINDING LOCAL BOUNDS 
ON THE SOLUTIONS OF BOUNDARY VALUE PROBLEMS, AND ON THEIR 
DERIVATIVES* 


By 
PHILIP COOPERMAN 


New York University 


1. Introduction. Methods for the approximate solution of the boundary value 
problems of mathematical physics are of little practical use unless there is reason to 
believe that the difference between the approximate and the exact solutions is sufficiently 
small for the purpose in mind. On the other hand, this difference, or error, is generally 
difficult to estimate. It is, therefore, a little surprising that there should exist a method, 
originated by E. Trefftz, which is capable of giving not merely an estimate, but upper 
and lower bounds on the solutions, and their derivatives, of the boundary value problems 
of physics. The method has been developed further by many writers. The purpose of 
the present paper is to encompass the previous work in a general scheme. In doing this, 
the writer borrows much from his predecessors; the debt is acknowledged by the bib- 
liography at the end of this paper. 

Before going into the general, and consequently abstract, theory of this method, we 
shall illustrate it by using it in the field of electrostatics. Let u, u, , uw, be functions of 
the rectangular coordinates x, y. Let D be a region, not necessarily simply-connected, in 
the z, y plane and let y, and y. be portions of the boundary of D such that y, + 72 is 
the entire boundary, y. Then, the general two-dimensional boundary value problem of 
electrostatics may be indicated by the equations: 


IU, > ‘ 
a + ous +f=0 in dD, (1) 
Ox Oy : 
OX 9 : 
U, oo Uy oy +gqg=0 on a (2) 
on on : 
u=k on Yo; (3) 
ou Ou F * : 
“u= ax’ U2 = ay in D+ 7’; (4) 


where ds is the differential of the arc-length along y. If u, u, , wu, denote the functions 
or Y2 as the entire boundary, thus eliminating equation (3) or (2), respectively. Also, 
the function f, which is 4 7 times the charge density, may be identically zero without 
invalidating the theory. 

The preceding problem, which may be labeled Problem I, is known to be equivalent 
to the following variational problem which may be called Problem I’. Let wu’, uj , uz be 
continuous functions satisfying equations (3) and (4). Then, Problem I’ consists of 


*Received January 21, 1952. 
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finding the particular set of functions of this type which minimize the energy integral 


given by 


. 


J,{u’] = Hi i} (ul? + us’) — ju} dx dy + | gu’ ds (5) 
D 


Y.1 


where ds is the differential of the arc-length along y. If u, u, , u2 denote the functions 
satisfying (1) to (4), then the statement that Problems I and I’ are equivalent amounts 


to saying that 


Min J,[u’] = J,[u]. (6) 


Moreover, by an application of Green’s theorem, we obtain 


: 1 f fe ; j Ox JY - 
J(u] = 9 es I fu dx dy + | gu ds + | Hw, = + Uy ou) ax (7) 


D 


Equation (6) is nothing else but the familiar Rayleigh-Ritz method, for it states 
that for any functions wu’, uj , ui of the class previously defined, the energy is greater 
than the energy due to the true potential. In order to obtain a lower bound for J,[u], 
we shall have recourse to a method similar to the Rayleigh-Ritz. Let us consider con- 
tinuous functions uj’, u3’ which satisfy equations (1) and (2). Note that therefore u;’ 
and u;’ are not necessarily derivatives of a function w’’. For this class of functions, let 
us consider the problem, which we shall call I’’, of extremizing the integral K,[u’’] de- 
fined by 


oar l é . . t Ox 0 4 
K,[u’’] = - | (ut’” + us’) dx dy + | kul’ — + ut’ oy ds. (8) 
2 Js ; : ‘ an “ an 
7 ve 
D 
K,[u”’] is obtained by forming the expression 2J,[u] — J,[u’] and by then substituting 


47 


for the functions u, uw’, etc., the functions u/{’, us’. 
With regard to Problem I’’, the following theorem can be proved. 
Theorem: 

(a) K,[u’’] attains a maximum value; 

(b) Max K,[u”’] = Min J,[w’]; 

(c) Max K,[u’’] is attained for the same functions u, u,; , u. which minimize J,[w’]. 
A discussion of the physical meaning of this theorem would be in order at this point. 

Let us consider, first, Problem I’. In this problem, the class of allowable trial functions 
consists of functions which satisfy the same boundary condition as the potential of a 
charge distribution given by the function f. In calculating the energy integral, J,[w’], 
one uses the derivatives of u’ as the components of the field strength. However, no 
attempt is made to have these derivatives satisfy the conditions imposed on the field 
strength (equations 1-2). Thus, the trial functions of Problem I’ satisfy the conditions 
imposed on the potential, but not those imposed on the field strength. 

In contrast to this, the trial functions allowed in Problem I” satisfy the conditions 
imposed on the field strength but not, except for the solution itself, those on the potential. 
Thus, in Problem I’, one attempts to approximate the potential whereas in Problem 
I’’, one attempts to approximate the field strength. 
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To prove the theorem, we start from the identity 


( ; 
r Jd 99 OF r, OF : 
| 1" (u % + us : 20) as “#% = (2e’ eet’) — +2 wut} dx dy. (9) 


By use of equations (1) to (4), this can be transformed into 


—| gu’ ds + | Kut! in ur 2H) ds = --|f fu’ dx dy 
on Or 


"7: ~*~ Ve 
D 


+ | (ufut’ + usus’) dx dy. (10) 


D 


But 
ulut’ + ubus’ = BA (ul? + ud?) + (ul? +0”) — (a — WP — (uh — t’)*}. (11) 
Substituting (11) in (10) and comparing the result with equations (5) and (8), one 


obtains 


es ° ’? l 4 , rr\2 f 2 
Jilu’] — Aylu”’] = 5 | iat — a) + Oe > uy’) } de dy. (12) 
D 
The right-hand side of (12) is never negative and vanishes if and only if 
du’ = : ou’ 
uy = — = uf’ ul, m= — = uy’. (13) 
Ox OY 
The rest of the proof is easy, for: 
(a) from (12), it follows that Min J,[u’] = A,[u’’], and hence K,[wu’’] is bounded above; 
(b) since Min J,[u’] = J,[u] and Max K,[u”] = K,[u] and since u, u, , U2 are allowable 
trial functions in both Problems I’ and I’, it follows that in this case (13) will be 
satisfied and consequently by (12), that J,[u] = K,[uJ; 
(c) this part follows from what was said above. 

This theorem enables us to bound the energy from above and below. In order to 
bound the potential function, therefore, we need to be able to express the value of the 
potential function at a given point as a function of certain energies. We start by con- 
sidering the expression for the potential involving Green’s function. This is 


er 


uzy=- I § f(é, n(o log r + r)} dt dn 


de [ {ols log r + r)} ds + [ {k 2 (2 log r + r)} ds (14) 


D 
where r = [(xz — £)* + (y — n)’]'” and R is a function to be discussed subsequently. 
Since f, g and k are known functions, we can evaluate o(z,y) where o is defined by 


I —— L Ff sa Rf, ‘i 
a(x, y) = = I S(E, n) log r dé dn + on a g log rds + ~ / k a (log r) ds. (15) 


D 
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Hence, 


, ° fr | OR ae 
u(z, y) — o(2, y) = — [| f(é, n)R dé dn + | gR ds + | k — ds. (16) 
JJ * ne day OM 
D 
Now, RF is the regular part of the Green’s function and so is a solution of the boundary 
value problem denoted by the equations 


0 ] P e 
—R,+—R, =0 in OD, (17) 
Ox OY 
Ox. OY l a "0 
> a s = zo yr = 8 
R, an + R, a -t- on an (log r) 0 on ">; (18) 
] 
R, = —=— logr on 2 (19) 
lls 
OR OR ; 
So ee Dee ee 
R, = ar? R, ay in D + y. (20) 


But this problem is of exactly the same type as Problem I and the corresponding varia- 
tional integrals, which we shall call J,[R’] and K,[R’’] may be obtained by setting f = 0, 


g = (1/2m) (0/dn) (log r), k = —1/2z log r in equations (5) and (8). Similarly, we have 
— 1ff a , , aR , \ 

Min J.[R’] = Max K.[R’’| = - , rr) ds — log r — ds?. (21 

in J2[R’] = Max K,[R”’| = 7- u. R + (log r) d in og 5, Oy (21) 

Let g = u — a R, where a is a parameter whose value remains to be determined. 


Then, g is a solution of the boundary value problem denoted by the equations 


< a+ “s gqa+f=0 in D, (22) 

q: - + ou +9-5 < (logr)=0 on 1, (23) 
q=k+t a log r on Ye > (24) 

“= of, = . in D+y¥. (25) 


Again, this problem is of the same type as Problem I and the corresponding variational 
integrals which we shall call J,[9’] and K;[q’’] may be found by making the substitutions 
g — a/2nx d/dn (log r) for g, k + a/2x log r for k. Hence, we have for J;[q] = K;[q], the 
expression 


= i | fq dx dy 


vd 


D 


Min J;[q’] = Max K;[9"’] = | 


> j eae ¢ z , 4 . a oq - 
¥ ] (0 2r dn log r)a “+ | (i T 2r log ?) on is. (26) 


- Fs Ys 
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Now, equation (25) can be put into the following form 


Min J;[q’] = Max K,[q’’] = 1 — [| fu dx dy 


D 


+ [ gu ds + I k= » 2u ash 


a si f , , 
+ ef ff FR dx dy — | gR ds — k a ash 
D 


ltt 2% lf aR } 
+ : it i R re (log r) ds — on , log r an ds 

a :. ¥ re) 1 f ou as 
+ 2 on a ua (log r) ds + on / i log r ax. (27) 


d, = Min J,[u’] = Max K,[u’’], 
d, = Min J,[R’] = Max K,[R’’], 
d,; = Min J;[q’] = Max K;|[q’’]. (28) 
Then, by use of equations (7), (16), and (21), equation (27) becomes 
f : ¥ 0 a i Ou } 
d, = d, +a’d, — = — , rr) ds - — rr— ‘ 9 
1, i tad, 9 \ o+ on I ua (log r) ds on I log 7 an ds (29) 


But, by Green’s identity, we can write 


1 f ee ee Oe - If {2 aR 4 ou oR | 
- in ua log r) d: on l= log rds = aw + by ay dx dy 
a u r) oR 
= If r+ +> Mt) a dy + I gR ds + [ k a ds 
fr | OR 
= - ff fR dx dy + [ gR ds +  & k i ds 
D 
=u-— do. (30) 
Hence, (29) reads 
d, = d, + a’d, — a(u — oa) (31) 
or 
atu — co) = d, + ad, — d;. (32) 


(32) gives us the value of the solution in terms of a known quantity o, a parameter a, 
and the minima, d, d, , d; of three variational problems. 

Since the quantities d, d, , ds; , are bounded above and below, we can also bound 
uw — o. Thus 

K,|u’’] + o’K.[R”] — Js[q’’] S a(u — 0) S J, [u’] + a°J2[R’] — Ks[q’’]. (33) 
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Furthermore, since q’ = wu’ — a R’, q”’ = u” — a R”, we can write 


Js[q’] = Jilu’] + aJ,[R’] — 2ab’, 


K,[q’’] = K,lu’’] + a’ K,[R’’] — 2ab’’, (34) 
where 
, l tj ‘Dr ‘pr vi “Dr 
b’ = 54 I| uiR{ + uzR) dx dy — I| fR’ dx dy 
J ] 0 ) ) ai 
+ gk’ — u’ — [log r]} ds, (35) 
I (ot 2r an llog 7] } 
ifr ; £ Ox oy 
"== I (wR + us’RS’) dx dy — log Aw |. 42’ —] ds 
2 \JJ —_— d an |, ‘an? an 
r Ox 2 ' 
+ | A(R? ; + Ry’ ou) ass. (36) 
J 4 on on ) 
Let a = J,{u’] — K,[u’’] and c = J.[R’] — K,[R”]. Then, the inequality (33) may be 
written 
—a — ac + 2ab’ S alu — oc) Sat ac + 2ab’’ (37) 


which can be broken up into 
a+al(u — oc) — 2b’] + ac = 0, (338) 
a+ al2b’’ — (u — o)|] + ac = 0. 39) 


Since (38) and (39) must hold for arbitrary values of a, it follows that their discrimi- 


nants must never be positive. Hence, 


tac — [(u — o) — 2b’]° = O, (40) 
4ac — (26 — (u— a)|° = v; (41) 
or 
2b’ — Vac) Su-—-o S 2b’ + Vac), (42) 
2b’ — Vac) Su-—o S 2b" + V ac). (43) 


Furthermore, both (42) and (43) must hold simultaneously. Thus, if b’ 2 6”, we have 


2(b’ — Vac) Su—o S 2b" + Vac) (44) 


and if b’ => b’, (44) holds with b’ and b” interchanged. In either case, we can define a 


“best’’ approximation to u which we shall call u* by 


u* =o + b’ + b” (45) 


and, it is clear that 


a= we? | 5 BVee — | 8’ — 6" | (46) 


(46) gives the bound on the “error” for the best approximation, 
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Examining equation (14), one notices that one can obtain a formula giving du/dxr 
and du/dy by straightforward differentiation. Since in this formula f, g, and k are functions 
of &, n and not of z, y, it is only r and R that is differentiated. Hence, one needs to set up 
variational problems for @R/dz and dR/dy, but not for du/dx or du/dy. Thus, Problems 
I’ and I” remain unchanged. Problems II’ and II” are obtained by differentiation of 
equations (17) to (20). Problems III’ and III” are set up in the same manner as pre- 
viously for a function g = u — a dR/dzx or gq = u — a OR/dy and, finally, one obtains 
equations entirely similar to (45) and (46). 

2. Notation. We shall now try to generalize this method so that upper and lower 
bounds on the solutions and their derivatives of a certain class of boundary value 
problems may be obtained. The class of boundary value problems which will be treated 
consist of all those having the same solution as a positive definite, quadratic variational 
problem. We shall here restrict ourselves to the boundary value problems of systems 
of one or more second order partial differential equations. However, the plate equation 
can also be treated by the same method. 

Let vy; and ¢;,(i = 1---m,k = 1 --- n) be functions of the independent variables 
x,, and let E’[y] represent a homogeneous quadratic form in the ¢, and ¢,, with coefficients 
which are functions of the z,. E’[¢] will be either positive definite or semi-definite in 
the algebraic sense, that is, E’[g] = 0 for arbitrary functions ¢; , ¢,, . The quantities 
E'{g] and E!,[{¢] are defined by 
dE’ |e] 
a ag 


aE'le) © 


E!,[¢] -_ Oy . 


Eile] = 


Using the equations above, two other expressions, L,[¢] and 1/,[¢] are defined by 


Lile] = Euale] — Eile], Mile) = Edlela.. (2) 
where 
” _ 9 py, _ 9x, 
Et, ele] —_ Or, Eile), Xe oo an’ (3) 


In these equations, n stands for the outward drawn normal to a domain D, and the 
summation convention is used. 

By E[¢], we mean the integral of E[g] over a domain D. Ely, ¥] and E’[¢, ¥] are the 
bilinear expressions corresponding to E[¢g] and E’[g]. We also define other integrals H and 


B by 
Hise = | --- | fea, (4) 


a oo / bes / go dA, (5) 


where y is the boundary of D and dV and dA are the appropriate volume and surface 
elements. B,[g, ¢] and B.[g, ¢] designate the integral Blg, ¢] extended over a portion 7, 
or y2 of the boundary, and the case y, = 0 or y2 = 0 is not excluded. In any case, we 
shall always have 7, + v2 = v¥. 

3. Some lemmas. In this section, we state without proof some simple lemmas 
which will be used frequently in the sequel. Although simple, these lemmas form the 


foundation on which the theory of the Trefftz method rests. 
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Lemma I. 
E'le, ¥) = ttEilelvi + Eilel¥s}- 
Lemma 2. Let ¥;., denote the derivative of y; with respect to x, . Then, 
E'ly, ¥] = 2{-¥iLile] + Wa — vi JEL] + WEale)) cS. 

Corollary. 
Ely, ¥] = 3{—A[L le], vi] + H[Eble], vie — Yin) + BIM Le], vil} - 

Lemma 3. The equations E’[y] = n; , E::{¢] = i, always possess solutions under the 
condition that if any pair of the E{{[¢] or E%,[¢] are identically equal, then only one 


of the corresponding equations should be used. 


Let ¢; , ¢;, be a set of continuous functions satisfying the conditions 


Lemma 4. tC, 


L{¢g]=0 in D, M [vt] =0 on 7, 


If H[E?,{¢], v]) = 0 for all such sets of functions, then v = 0. 


4. The complementary principle. Let a set of functions u; , u,;; exist such that they 
satisfy the conditions: 

L;{u] = f; in D, (1) 
M lu] = —g; on > (2) 
ve; = £; on Y2 5 (3) 

Ou; ‘ a 
“ij = Ui =o in D+ y. (4) 

Od 


The problem of finding these functions u,; , u,; is the fundamental boundary value 
problem which we shall try to solve. We shall call this Problem I. 
Let V’ be the class of all continuous* functions u/ , u/; satisfying equations (3) and 


(4). Let the functional J,[w’] be defined on the class U’ by 
Ji{u’] = Elu’]) + Hf; , ut] + Bilgs , ui). (5) 


The problem of finding those functions u’ , u’; for which J,[u’] attains its minimum is a 
7 
variational problem which we shall call Problem I’. By the Dirichlet principle, we know 
that Problem I and Problem I’ are equivalent, i.e., they have the same solution. 
Let U”’ be the class of all continuous functions u?’ , ui; satisfying equations (1) and 
2), and let the functional K,[w’’] be defined on U” by 
1 J 


K,[u’’] = —E[u’’] + B,[M,[u’’], ki). (6) 


The problem of finding those functions u‘’ , u/} for which K,[u’’] attains a stationary 
value is a variational problem which we shall call Problem 1”. 

The complementary principle is expressed by the following theorem. 

Theorem I. 
(a) Problem I” is equivalent to Problems I and I’. 


*For the sake of simplicity, we shall consider only continuous functions. However, less restrictive 


conditions could be imposed, e.g., piecewise continuity. 


1953] SOLUTIONS OF BOUNDARY VALUE PROBLEMS 367 


(b) K,[u”] attains a maximum value. 
(c) Max K,[u”’] = Min J,[u’] 

Proof. (a) Let v; and v,, be functions of the class U” and let ¢; and ¢,, be a variation. 
Thus, the ¢’s satisfy the equations 


L.{t]=0 in D, M,[t]=0 on y,. (7) 


If v; and v,, are to provide a stationary value for K,[{w’’], it is necessary that they satisfy 
the condition 


—2E|v, ¢] + B.[AT,[¢], ki] = 0 (8) 
for all ¢. By the corollary to Lemma 2, this condition becomes 
AEGAS), v6.4 — vue) + BAM [E], ki — v4] = 0. (9) 
It is clear that we can choose ¢’s such that M,[¢] = 0 ony... Hence, by Lemma 4, applied 
to (9), the conclusion is 
V5 = Vn (10) 
Equation (9) now reduces to 
B,[M [fg], ki — v1] = 0. (11) 


It can be shown that ¢’s can be found such that the /,[t{] are equal to arbitrary con- 
tinuous functions. Hence, 


v, = k; on Yo + (12) 


Thus, the v, and v,, satisfy all the conditions of Problem I. This proves Part (a) of the 
theorem. 

(b) Since the v; and v,, are solutions of Problem I, and since we assume that solu- 
tions of Problem I exist, we know that the »; and v,, exist, and we can identify them 
with the u; and u,;, . Hence, K,[u’’] actually attains a stationary value. Because of the 
fact that E’[w’’] is never negative, this stationary value must be a maximum. 

(c) By the corollary to Lemma 2, 

Efu] = 4{-H[f; , uc] + Bilg: , uc] + B.[M,[u], k,]}. (13) 


Then, 


K,[u] = Max K,[u’”’] = 3{H[f, , uw) + Bilg: , us] + Bo[M fu], k}}, (14) 


and also 
J,{u) = Min J,[w’] = 2{H1f; , ue) + Bilgs , us] + BLM. [u], ki}. (15) 
Comparing (14) and (15), we have 
Max K,[u’’] = Min J,[u’]. (16) 


This completes the proof of the theorem. 
5. The auxiliary problems. Let F;; denote the components of the fundamental part 
of Green’s matrix and R;; , the components of the regular part. It is assumed that the 
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F,; are known functions. Then, for a given value of j, say 7 = k, the functions R;, solve 
the boundary value problem indicated by the following equations: 


L;([R,] = 0 in D, (1) 
M{R.]J= —-M[Fi] on mn, (2) 
Ry = —Fy on Yo; (3) 
Ru; = Ri; in D+y¥. (4) 


This is a problem of exactly the same type as was treated in the preceding sections. Let 
us call it Problem II. 
Problem II’ consists of finding the functions which minimize the functional J.[R/] 


defined by 
J[Ri] = E(Ri] + B,[M,[F,], Ri] (5) 


among all continuous functions R/, , R?,; which satisfy conditions (3) and (4). Problem 
II’ consists of finding the functions which maximize the functiona] K,[Rj’] defined by 


K,[Ri’] = —E[R;’] — B[M.[Ri’], Fit] (6) 


among all continuous functions R’/,;, R’/,; satisfying conditions (1) and (2). 
In the same way, we can consider the problems which have as their solution the 


functions, Q;, = u; — a Ry , where a is a parameter. The boundary value problem 
which we label Problem III is indicated by the equations 

L;{[Q.] = f; in D, (7) 

M [Q.] = —g: + oM,[F;,] on v1; (8) 

Q. =k + aF,, on 2, (9) 

Qa; = Q.; in D+ y. (10) 


Problem III’ consists of finding the functions which minimize the functional J;[Q;] 
defined by 
J3(Qi] = E(Qi] + Afi , Qe) + Bilgs — aM [Fi], Qe] (11) 
among all continuous functions satisfying (9) and (10). Problem III” consists of finding 
those functions which maximize the functional K,[Q;’] defined by 
K,[Qi’] = —E([Qi’] + B.[M.[Qi’], ki + oF iu] (12) 
among all continuous functions satisfying (7) and (8). 


6. Bounds on the solution. If u, denotes the value of the function u,(z,y) at the pole 
of the Green’s matrix, then it is well known that 


“oO, = Hif; ’ R| + Bg; ’ Rx] i B,([M [Rx], k,), (1) 
where 
o. = Alf; , Ful + Bilg: , Paul + B.[M.IF;], k;). (2) 


Our object at this point is to express u, — o, in terms of the stationary values of the 
variational problems set up in the two preceding sections. 
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Consider the stationary values achieved in Problems III’ and III”. By the comple- 
mentary principle, these stationary values are equal. Since Q;,, = u; — a R,, , we can 


write 
J;{Q.] = K3[Q,] = Elu] — 2ak[u, R,] + a°E[R,] 

+ Alf; ,u:) — eAlf;, Rul + Bilgi , us) 

— aB,|M,{F,], us] — aBy lg: , Rian) + a Bs[M[F,), Rie] (3) 
which by considering equations (4.5) and (5.5) can be transformed into 


J IQ, | J, {u] - a’ J.[R,] _ a{2E[u, R,] + H{f; > Rix] 


+ B,[M,[F,], ui] + Bilgs , Ril}. (4) 
3y the corollary to Lemma 2, 2 E[u,R,] can be transformed. This gives, taking ¢; = 
RW = U; 

2E[u, R,] = B[M,[R,], us] = —B,[.[F,], us] + B.[M,[R,], ki). (5) 
Consequently, we have 
2F lu, R.) + Alf; , Ral + BMAF), us] + Balg: , Rul 
= Hf, , Ral + Bilgs , Ric] + Bo[M.(Ri], ki) (6) 


which by equation (1) is equal to u, — o, . Hence, (4) becomes 


a(u, — ox) = J,[u] + a’J2[Ri] — JslQx]- (7) 
Let 
a = J,{u’] — K,[u’’] = 0, 
bi = J[u’, Ri], 
bi’ = K,[uj’, Ri’), (8) 
c, = J,[Ri] — K,[Ri’] = 0, 
where 


2J,[u’, Ri] = 2E[u’, Ri] + Hf; , Ri] + BAF: ), uw) + Bilge, Rial, (9) 


and 
2K,[u’’, Ri’) = —2E[u’’, Ri’) — B.[M,[u’’], Fu] + B.[M [Rv], k,). (10) 
Then, since each of the quantities on the right-hand side of equation (7) have upper 
and lower bounds, e.g., K,[u’’] < J,[u] < J,[w’] (7) can be turned into the following 
two inequalities. 
a+ {(u.— o) — 2bi}a + ca° = 0, (11) 
a a {2b;’ zoe (uy = a.) }a “bh c,a° = 0. (12) 
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Since the discriminants of the quadratic forms on the left of these inequalities must be 
non-positive, we arrive at two further inequalities. 


2(bi — Vac,) Su — o S 2bi + Vac): (13) 

2(bt’ — Vac) <u — ox S 2(bf’ + Vac). (14) 
Since both inequalities must hold simultaneously, we have, if bi S by’ , 

2(bi’ — Vac) S um — 1 S 2Abi + Vac), (15) 


and if bj = b;’ , the same inequality with bj and b;’ interchanged. 
It follows from (15) that the best approximation, u% to uw, is given by 


ut = o, + BL + Df’. (16) 


For this choice of u*, we have that the “error’’, namely, | u, — u% | , is bounded by 


lu —ut| Ss 2V ac, — | bf — bf’ |. (17) 


It is interesting to see the meaning of the quantities, bj and bj’ , on the right-hand 
side of (17). By applying the corollary to Lemma 2 to equation (9), we get 


2b) = 2J,(u’, Ri] = Alf; , RE) + Bilg: , Ri] 
+ B,[M,{RU, k.] — HILARY, uw!) 
+ B,[M,[F;] + M,[Ri, uf]. (18) 


The first three terms represent an approximation to u, — o, in terms of the functions 
R’,,. The last two terms would vanish if the PR’, were actually the regular part of the 
Green’s matrix. Hence, 26; is an approximation to u, — o, . A similar interpretation 
can be given for 2b,’ . Thus, u% is actually the sum of an “exact term’’, o, , and the 
average of two approximations to the remainder, one from above and one from below. 

7. Bounds on the derivatives. By differentiation of equation (6.1) with respect to 
the zx; coordinate of the pole, we get the relationship 


Uu.; — or; = Alf:, Rui] + Bilg: , Ru.j] + Bo[M[R..,], &i). (1) 


The boundary value problem solved by the functions R,,,; (&,7 fixed indices) can be 
derived from Problem II by differentiation of equations (5.1) to (5.4). This gives 


L,{R:,;] = 0 in D, (2) 
M.IR;.;)= —MiAF.:;] on 7, (3) 
Ru = —Fis.; on V2» (4) 
Riss = Riza; _ in D+ y¥. (5) 


Thus, the functions R,,.; satisfy conditions analogous to the conditions of Problem II. 
Then, just as we did previously, we can define upper and lower variational problems. 
The same can be done for the functions Q;,.; = u; — a Ry,; . The theory carries through 
in exactly the same way and enables us to find a best approximation u%,; and bounds on 
the quantity | u%,; — u%,; | 

8. The boundary value problems of elasticity. As an illustration of the application 
of this method, we consider the equations of elasticity. We identify the u, as the dis- 
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placements. If e;; represents a component of strain and S,; a component of stress, then 
these quantities are defined by 


Qs; = Fuss + Uj), (1) 
Si; = CijeiCt » (2) 


where the c;;,; are symmetric. The boundary value problem is given by the equations 


L,{u] = Si, = f. in D, (3) 

M [ul] = Si;25.. = 9: on Y ; (4) 
u; = k; on 2, (5) 

Ug = Us; in D+y. (6) 


The quantity E[u] is defined by 


ial = alll €,;8.; dV. (7) 


The quantities a, bj , bi’ , c, are given by 


= 5 fff esi. av + Ii fut dV + If err 
+3 /f/ eft dV - | Siiz,.k:dA (8) 
2b, =} fff ter,S.sthe] + eulRASt} av : 
; 
+ II fin dV + I gi, dA + I S.[FsJa;..uidA, (9) 
. " os 
ol “all ef) S,,[Ri’] + e:;[Ri’|Sii} dV 


+ If» S.[R']t;..dA — [f Fast .ndA, (10) 


CG = 2 ARES, ;[Ri]} dV + Si;[F iz; nRi, dA 
a If I 


+3 fff fe,,{Ri)S.s(RV1} dV + il] Si(RltinFu dA, (11) 
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where e;; , S{; mean that these quantities are defined by equations (1) and (2) replacing 
u;; by uj; . Expressions such as e;,,;[Ri], etc., have a similar meaning. It then follows 


that the best approximation to u, is u% where 





ut 


and, the error is bounded by 


The stresses 
derivatives of the u; . 


9. Boundary value problems of single, second order, elliptic equations. 


(12) 


o, + bi + by’ 


2V ac. — 


(13) 


IIA 


bi — bi’ 


and strains can then be bounded by applying the same method to the 


It is well 


known that any linear, self-adjoint, elliptic partial differential equation can be put in the 


form 
Au — du = f (1) 
Uniqueness of the solution is guaranteed if we take d(z,y) > 0. The boundary value 
problems associated with this equation can be written 
ra] 0 , E : 
L ie a ee = Oy ae 7 in res (2) 
or OU . ; 
Wlu] = ua, + uy —g on 7 (3) 
u= } on Y ( t) 
Ou Ou ‘ 
l; ee ilo = in D + Y (5) 
Ox oy 
E[u] is defined by 
Blu > | (ua; + us + dw) dx dy. (6 
D 
In this case, the qua ies a, b’, b’’, c are given by 
Lore rp 
a 5 || - / + ax dy — | Ju dx dy 
f a l rf 472 pp? roe i a ‘ a rT \ _ 
+ | ( as + > (U T Uz + au ea. te — | INCH Zz. a U2 Yn) ds, (4 
2b’ = |] (WR! + wR} + du’ R’) dx dy 
dy + | gh’ ds + | (Fiz, + F.y,)w’ ds, (8) 


I fR’ dx 


D 


J ¥ 
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20" = —|f (witht! + ufRY + dw’ R”) de dy 
"D 
+] KRYr, + Ry) ds— | Flute + uly) bs, 
“es * Te 


¢= z I| (R? + RY + dR”) dx dy 


D 


+L [ (Fix, + F.y,)R’ ds 


+5 ff (ae? + Re? + aR’) de dy + | 


ma” 


(Ri'xa, + Ro’y,)F ds. (10) 


D 


Then, just as before, we have 


ut=ot+b’+ bd”, (11) 
u—u*| < 2Vac — | b’ — b” (12) 
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ON MAPPING AND MEASUREMENT OF RANDOM FIELDS* 


By 
MAHINDER 8S. UBEROI anv LESLIE 8S. G. KOVASZNAY 
The Johns Hopkins University 


1. Introduction. [Tor random fields in space, such as turbulent flow patterns, the 
significant measurements are presently believed to be those of such statistical properties 
as the correlation functions and the “‘power spectra”. Such measurements are usually 
made difficult by the fact that we can seldom measure directly any quantity at a point 
where, by point, we mean a region of space so small that the properties are substantially 
constant within the region. When the resolving power of the instrument is comparable 
with the extent of the fluctuations occurring in the field, a correction must be made for 
the resulting averaging that will inevitably occur. A further complication frequently 
occurs in that we are not measuring the desired quantity itself but some functionally 
related quantity involving derivatives, integrals and the like. Again it is necessary to 
correct the statistical results for such unwanted effects. By a proper choice and use of 
instruments it is usually possible to obtain a linear response from the instrument. For 
a very large class of useful instruments, the instrument reading Q(x) is expressible in 
the form 


Q(x) = | K(x, s)U(s) dV(s) 


where U(x) is the quantity whose measurement is desired and the kernel A(x,s) is 
determined by the properties of the instrument. x and s are space position vectors and 
dV(s) is the volume element at the point s. Unless otherwise noted the integration is 
extended over the entire field. A few simple examples will perhaps prove useful in illus- 
trating the generality of this type of response behavior. 

1. If the instrument is perfect, then K(x,s) = 6(s — x) where 6(s — x) is the Dirac 
function. In this case Q(x) = U(x). 

2. If the instrument measures a derivative of U(x), then A(x,s) = —6’(s — x) where 
5’ is the corresponding derivative of the Dirac function. 

3. If the instrument gives a “Gaussian average” over a region of space having a 
characteristic dimension a, then 


K(x, s) = —n exp [=== 8] 


a 


where n is the number of dimensions of x. 

We can view the above integral equation as a mapping process in which the instru- 
ment performs a linear operation, converting the original field U(x) into the associated 
field Q(x). From the instrument readings, we have the map of 2(x), whereas the map 
of U(x) is the one that existed in nature. For random fields, the problem is to determine 
the statistical properties of U(x) from the statistical properties of 2(x). 

In what follows, we have been able to solve this problem for the case where the field 


*Received February 25, 1952. 
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U(x) is statistically homogeneous infinite in extent and the response characteristics of 
the instrument are not changed when the instrument is moved around in the field. 
Because of this independence of position, the response kernel must be a function of 
the single variable (s — x). Thus we have 


oz) = | K(s — HU(s) aV(s). 
A simple change of variable gives the equivalent form,* 
ox) = | K(s)U(x + s) dV(s). 


2. General relations. For definiteness, and with application to turbulence in view, 
we consider a three-dimensional random vector field, U,(x), which is statistically homo- 
geneous and infinite in extent. 2;(x) is an associated random vector field which is con- 
nected with the primary field by the relation 


0.(x) = | K,,(s)U,(x + s) dV(s) 


d 


where a repeated index means summation over the index. Averages are taken over the 
entire range of independent variables. 


(Us@) ue. = lim > | Ue) aV@) 


aan, VE) ity 


where V(r) is the volume of a sphere of radius r. We assume (U,(x)),v. are all zero. 
Space averages may be functions of a parameter like time, however, we omit reference 
to it to avoid an elaborate notation. R,;(&) and 8,;(&) denote correlations of the origina} 
and the mapped field respectively. 


Ri(E) = (U(x) U(X + B))av. ; 
B:(E) = (2;(x)2;(X + E))ev. ; 


E,;(k) and T,,(k) are the spectra of U,(x) and 2,(x) respectively. 


E,,(k) = = | R,,(&) exp [ck-&] dV(é), 
Sr 

r,,(k) =< J / 8;,(€) exp [ck-E] dV(é), 
8x 


with corresponding inverse relations, 


R,(®) = fz. (k) exp [—7k-E] dV(k), 


B; ;(E) [ I’, ;(k) exp [—7k-E] dV(k). 


*These operators have been used, among others, by N. Wiener). 
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It is well known that the auto-correlation (e.g. R,,(&)) and the “power spectrum” are 
the Fourier transform of each other (1). By power spectrum of U,(x) we mean the 
quantity, 
|2 
mer © [ — U,@ exp [ak-x] aV(x) | 
“Vi(r) | 
where V(r) is the volume of a sphere of radius r. This limit is identical with £,,(k) which 
is the Fourier transform of R,,(&). 2,,(k), E22(k) and E;;(k) are the power spectra of 
the vector random field U,(x), similarly I',,(k), T'..(k) and T;,(k) are the power spectra 
of 2,(x). In general, by power spectra we mean the diagonal elements of F,,(k) or I, ;(k) 
and by spectra we mean the entire tensor. 

For simplicity we have assumed pure band spectrum. In order to include the case 
of line spectrum, it is necessary to define an integral spectrum. The correlation is the 
Fourier-Stieltjes transform of this integral spectrum. 

As mentioned earlier, 2;(x) can be viewed as the mapping given by a linear instrument 
which is exploring the field, U;(x). Since Q,(x) is the field directly available to us from 
the output indicators of our instruments, the usual situation is that 8;;(&) and T,,;(k) 
will be known to us by means of the conventional experimental techniques for obtaining 
such statistical data. The problem before us is to obtain the statistical quantities R, ;(&) 
and E;;(k) of the primary field if we know the response function K,;(s) of the instrument. 

If we substitute the value for Q; = f K,,(s) U,(x + s) dV(s) into the expression 
8.;(—) = (Q;(x) Q,(x + &)),.. , we obtain 


B.,() = lim vA | avi) I| K,(t)K;(s)U(x + U(x + & + s) dV(t) dV(s). 


) . 
JV(r 


The integration involving an averaging process and the two integrations with respect 
to t and s can be interchanged if K;,;(s) is zero outside of a finite domain or if U,(x) is 
finite 


a.) = |f Kuk (= Lim ro |, Ue + OUI TEF ave | dV(s) aV(t) 


= I Ki. (OK; (s)Ril— + s — t) dV(s) dV(t) 


let 


then 


B;,(§) = I Ku(s + 2)K,.(s)Rul€ — 2) dV(s) dV(s). 


Integrate first with respect to s and denote by Yix;:(t) the “correlation” of K,.(s) and 
K;,(s). 


a” [ Kuls + *)K,,(s) dV(s), (2.1) 
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then 
6;,(&) = [ Ving (DRE — 2) dV(e). (2.2) 


It is well known that the Fourier transform of the convolution of two functions is the 
product of their respective transforms. Taking Fourier transform of (2.2) we get 


Pik) = Sap) En(k) (2.3) 
where 
Siej(k) = / bck) mm [tea dV tO. (2.4) 


We can express E;;(k) in terms of I[;;(k) by solving the set of nine linear equations 
(2.4), and taking the Fourier transform of E,;(k) (which is now expressed in terms of 
T,;(kK)) we can get the correlation R,,;(&). So that equations (2.1)-(2.4) give the solution 
to our problem of expressing the correlations and spectra of the original field in terms 
of the corresponding quantities of the mapped field. 

If the matrix K;,;(s) is diagonal there is no cross feeding from one component to 
another. Each component can be treated independently and K,,(s) in effect becomes a 
scalar, e.g. for the first component we have 


7 


2,(x) = | K,,(s)U,(x + s) dV(s) 


and the formulas (2.1)-(2.4) simplfy to 


$40 = [ KuoKu(s + 2) dV(s), (2.18) 

B,,(é) = / Wile)R,,(E — 7) dV(«), (2.2a) 

r,,(k) = S,(®E,,(), (2 2a) 
and 

S.(k) = / di(2) exp [—ik-2] dV(2). 2.4) 
Now, 


¥:(t) = / K,,(s)Ki,(s + +) dV(s). 
A simple change of variable gives 
w(t) = | Ku(-s)Kule — 8) dV(s), 
ie. ¥,() is the convolution of K,,(— s) and K,,(s). If F(k) is the Fourier transform of 
K,,(s), then S,(k), the Fourier transform of ¥,(+), is F( —k) F(k) = | F(k) |’. Therefore 


S,(k) is a real and an even function of k. We may call it the power sensitivity spectrum 
of K(s); it gives the square of the sensitivity of the probe as a function of the vector 
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wave number, k. It is a multi-dimensional spatial generalization of the customary 
frequency response in circuit analysis. The principal difference, is in the fact that time 
has a preferred direction, on the other hand space does not necessarily have such a 
restriction. In the standard methods of circuit analysis K(s) is the impulse response of 
the circuit and it acts only on the “past’’ of the signal. 

Before taking up some cases in detail we consider the simple case of an instrument 
which has uniform sensitivity in a spherical domain. 

In the electron circuit theory the filter that averages the input signal over a finite 
interval in time with uniform weighing has no great significance. On the other hand, 
when instruments with finite spatial resolution are considered this is one of the simplest 
cases. A hot-wire anemometer with uniform sensitivity over a finite span is the one- 
dimensional case. A circular disc averaging some optical effect is a two-dimensional 
example. A small spherical probe measuring chemical concentration in“ fluctuating 
field can be a three-dimensional case. Fig. 1 shows the response function K(s), its ‘‘auto- 
correlation”, ¥(*), and its “power sensitivity spectrum’, S(k). Since the dimensional 
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form of these is different in one, two and three dimensions, these quantities are normal- 
ized so that f K(s) dV(s) = 1. 

3. Example 1: Velocity field measurements by hot-wire anemometer. We apply 
the considerations of the last section to isotropic turbulence. For this purpose we intro- 
duce the customary terminology. U;(x) is the fluctuating velocity field, R;;(&) its correla- 
tion tensor, and £,;(k) its spectral tensor. If the turbulence is isotropic R,;(&) can be 
expressed in terms of two scalar functions f(r) and g(r) (2) 


9 


(U; avg(7) = R,,(0, T, Q) 


(Utah (7) = Ri, (r, 0, 0) 


wt 


R,(&) = (U3 sO 5 ©) t+ o(r) 5: r” = Ef, . (3.1) 
lor J 
R;,(é) is an even function of &, and it follows that its Fourier transform E,,(k) is also 
an even and real function of k. f(r) and g(r) are connected by the continuity equation of 
incompressible flow (2). 
r of 
2 or° 


f(r) —- gr) = 


Total turbulence intensity 


> (U2. = R;(0) = / E,,(k) dV(k). 


E,,(k) is thus the spectral density of pe (U;)v. in the wave-number space. It is useful 
to introduce a ‘‘three-dimensional’”’ energy spectrum which will give the energy density 
of all disturbances having the same wave number magnitude k. This is done by inte- 
grating E,,(k) over a sphere of radius k, thus, 


E(k) = | E,:(le) do(k) 
where de(k) is the surface element of a sphere of radius k. If we apply the conditions of 
isotropy and continuity equation of incompressible flow to E,;(k) and make use of the 
definition of E(k) we find (3) 
E(k) 


 / — dé 2 —- 7 * 29 
B,(le) = ors (K°8is — kik). (3.2) 


Commonly measured “one-dimensional” spectrum is the Fourier transform of (U})... 
f(r). 


rep 


| if E,,(k) exp [—ik,r] dk, dk, dks 


(Ui)av.S(r) = R,,(r, 0, 0) 


= [ E.G) exp [—ik,r] dk, 


where 
1 ff B®) 


an ya (ke + ks) dk, dks 


E,(k,) = | | E,,(k) dk, dks = 

‘ (3.3) 
1 f° E(k) 
4J,, & ( 


a 


k* — k;) dk; k? = kk; 
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From which the inverse relation follows: 


) = 2x-| @ (1 a) | 
E(k) = 2k E k, dk, bok 


For the measurement of U,(x) component of the turbulent velocity, a hot-wire of 
length 2 1 is set parallel to z; . It is sensitive to U,(x) and its output is proportional to 
the integral of U,(x) over a length 2 / along the wire. For simplicity we have assumed 
uniform distribution of temperature along its length, giving constant local sensitivity 
to U,(x). 

The point of view adopted here is that simultaneous measurements of three-di- 
mensional random field U,(x) are made at all points in space and we represent these 
measurements by another random field 2, (x). It is from Q,(x) that we compute various 
correlations and spectra, and this step involves no “error”. The “error” is involved in 
getting 2,(x) and consideration of this “error” will give a “correction” not only for the 
correlations but also for the various spectra. 

In actual practice the wire is stationary with respect to the mean flow carrying the 
homogeneous turbulence and it is sensitive chiefly to U, up to a reasonably large tur- 
bulence level. Under the same assumption we can assume that a “frozen” turbulence 
pattern is swept over the wire and, knowing the mean speed, we can convert time fluctua- 
tions into space fluctuations. The wire resolution is effectively perfect in z, and 2, direc- 
tions, but is finite in x, . From this it is immediately obvious that it will map an isotropic 
field into a non-isotropic one, with the highest k; wave-numbers washed out by the 
integration (or “averaging’’) over 2 / at each wire position. In fact, only a probe with 
spherically symmetric sensitivity distribution transmits isotropy. 
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Figure 2 — Functions choracterizing the response of a hot wire of length 21 
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The mapping operator corresponding to hot-wire of length 2 / set parallel to 2; is, 


( &(s,) 8(ss) for |s|<J, 


K(s) = ) 
0 for | 83 | > U. 


We first calculate y(<) and S(k), the functions characterising the hot-wire or the mapping 
operator K(s) 


v(x) = | K(s)K(s + =) dV(s) 


ds; I 5(s,) 5(s2) 6(s; + 71) 5(S2 + 72) ds, ds, 


= 6(7,) 6(72)(2l — | 73 |); v% | XSI. (3.4) 
Since 6(s,) and 6(s,) are not functions in the ordinary sense, but infinitely sharp dis- 


tributions, we have shown the graph of ff K(s) ds,ds. and ff ¥(*) dr,dr, in Fig. 2. We 
see that while {{ K(s) ds, ds. extends over a distance 2 l, [f ¥(*) dr, dr. extends over a 
distance 4 J. S(k), the power sensitivity spectrum is the Fourier transform of (+) 


a2l a 
S(k) = | dr; ! 5(7,) 6(72) exp [—ik-*] dr, dr. . 


The graph of S(k) is shown in Fig. 2. It is constant with respect to k, and k, and it is 
almost zero for | k; | > 2/1. The area gibi the tails is not large and as] > 0 S(k) ~ 47 
for larger and larger value of k; and as 1 © S(k) acts almost like the Dirac function 
in x; direction. In these limiting cases it is necessary to properly normalize S(k). 

The relation between I’,,(k) and E,,(k) is given by (2.3) 


sin k, 
s(n )E Bulk) 


od (amt fy i) (2 + 


2r k 3 


I’, ,(k) 


for isotropic turbulence (see (3.2) ). I',(k,) is the ‘‘one-dimensional” spectrum of the 
mapped field. It is the apc, obtained by analysing with a wave analyser the output 
of a hot-wire sensitive to U, 


r,(k,) = il T’,, (kt) dks di 








(3 5a) 
saaes ~ || (in al)’ HK UK) (ke + k2) dk» dk, 
let c = (k; + k3)* and k, be the new variables. 
. 2° , kel ra 
T’,(k:) = - / do | i ” dks (ae yo CO agi 
(3.5b) 


ef wa? — 0/7] 2D a - Kak 
Oks 


II 
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where 
7» _ 2 [* (sin uy — 
W(a) = =| ( y ) @- yy (3.6) 


W(a) has been numerically integrated and its graph is drawn in Fig. 3; the dashed line 
is the asymptote. Comparing (3.5b) and (3.3) we see that for small / the measured “‘one- 
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Figure 3 — The function W(x) 


dimensional” spectrum IT, (k,) approaches the true “one-dimensional” spectrum EF, (k,) 
except for a factor 4 I? which is the steady state hot-wire response. For large I (i.e., 1 > ©) 
Tih) w [ A) (2 — 2)? dk. (3.7) 
l Jn, &K 
If the integral (3.5b) converges, then it is always possible to find an / such that (3.7) 
is correct within a preassigned approximation. A better estimate on the error in (3.7) 
will depend on k, and E(k,) for k > k,. The squared output of the probe is increasing 
continuously with increasing length of the hot-wire, for the limiting case of infinite wire 
length the squared output per unit length, tends to a limit. In order to put this in a 
standard form for inversion, we differentiate, T',(k,) with respect to k, thus: 


MiG) _ t E(k) __ dk 
lk, — i. ke (k? vm 9 aie 
For this limiting case the inversion of the integral equation is (see Appendix): 
, ow Zhi” por —_ ae __ 
EQ, = | Grae. (3.8) 


Integrating by parts, we get 


; 2k or? rk) — T'(k, 
E(W aos, & 2 (k) = I"(ky) 


xl (ke? — ky? k dk ° (3 8a) 
vk, \ 1 
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Integrating again by parts, we have 





2k; [ r’'(k, (kh, — k) — (Th) — TH) (2k? _ k?) dk. (3.8b) 


ek | 
( k=k, rl Jas (k? os ki)” 
Also 
ice 1 EW) 
E,(k) = if 4 (k® — Ki) dk 
—— a (k? — ki) Ts’) 
= Onl ; ah [ ds Ke _ sa 


integrating first with respect to k then integrating twice by parts, we get 


E(k) = 2 


_ a ; 


l [oe ‘(€) dé (3.9) 


We thus see that either the “‘three-dimensional”’ spectrum or the true ‘‘one-dimensional”’ 
spectrum can be recovered from the spectrum measured by a hot-wire of “‘infinite”’ 
length, (by “infinite” length we mean a length much larger than any scale of turbulence). 
For this limiting case the length to diameter ratio of the hot-wire is large and the tempera- 
ture distribution becomes strictly uniform and independent of the operating conditions, 
as we have assumed here. If 


E(k) ~ k* exp [—ki/ki], 
then 


27,2 
E,(k,) ~ exp [—ki/ks] and f() = exp | r 


where ky is a constant. In this case there is no relative distortion of the ‘‘one-dimensional”’ 
spectrum or the longitudinal correlation measured with a hot-wire of finite length. 





Calculations show that 


erf (kol) < 1 


4PE,(k,) lk, 
=~] for kol «1. 


T’,(k) Vr 


The relation between true and measured intensities of turbulence can be put in terms of 
the scale of turbulence (4,5). 


measure d (Ui) se - r _ i (z!) 
true (U?). a ef (hl) = a1 AT. 


where 


a iz f(r) ar. 


In these formulas we have taken account of the fact that the steady state response of 
the hot-wire is 4 1’. 
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The measured lateral correlation 8,,(0,7,0) is related to R,,(0,r,0) by (2.2a) 
bu) = [ ¥@RAG— 2) aV(e). 
Substitute the value of ¥(), from (3.4) 


= [ary ff a7) (7)(21 = | 2 DRE - 2) drs dry 


Str 
to 


Bis(E; ’ 


We consider two special cases of 8,,(&), 8,,(0,7,0) and 6,,(r,0,0). Substitute the value of 
R,,(& — +) from (3.4) we get 


81(0, 7,0) = Uae f(t | rs Dal? + + 73)'"] drs (3.10) 


and 


B,,(r, 0, 0) 





p2t 2\1/2] __ 1/2 
- (U?),, (21 — | - yale + a Be a — gl? + + 73) J 2 + gl(r? tp ayy} drs 
J-21 r + 73 
using the relation f(r) — g(r) = —r/2 f’(r) we get 
2l 
Bu(r, 0,0) = (Ui 2t ff? + 23)") drs « (3.11) 
Jo 


Equation (3.10) was first derived by Skramstad (6) by considering the relative 
geometry of two wires used to measure g(r). The main advantage of our point of view 
is that it is not necessary to consider separate corrections for spectra and correlations 
if we concentrate on the properties of the probe (linear operator). If the wire length is 
much larger than scale of turbulence, it may be considered effectively infinite, then 


B10, 7,0) pron aie 1/2 
WO = (U4 | gle? + 1)"] dr. 
y “0 


This integral equation can be solved and the true correlation recovered from the measured 


correlation (see Appendix) 


ae 1d : Y 8,,[0, (r? + iz 0) + 219 
U)w0) = Fa," | P+e " sadens 
let £ = r tan 6, then 
—i@ oie en 
(Un)w9(t) = Onl dr B,,(0, rsec 6, 0) dé. (3.13) 


4. Example 2: Density field measurements by shadow method. The shadow picture 
of projectiles at high velocities show vividly the turbulent structure of the wake (Fig. 
4); the intensity of light on the photographic plate is a function of the density of the 
medium surrounding the bullet. We want to investigate whether the ability of the 
shadowgraph to show details can be put to advantage in recovering the statistical in- 
formation about the density field from such a “picture” of turbulence. Firstly, the 
shadow method is sensitive only to the second derivatives of the density field and 
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Figure 4—The Shadowgraph of a 19mm. Shell at 1700 Ft/Sec. in Atmospheric Air. 


secondly, since the light passes through a finite thickness of the wake, the final change 
in the intensity of light is the integral of the effects along its path. So that the “picture”’ 
is only indirectly related to the turbulent density fluctuations. According to our termi- 
nology this picture is a mapping of the primary density field. We will first consider the 
relation between the light intensity on the photographic plate and the density field. 
The wake of a bullet is given as an illustrative example and we consider below the 
random density field in general. 

Let p(x) be the fluctuating random density field, statistically homogeneous and 
infinite in extent. A portion of the field, lying between x, — 1 and 2, + / and extending 
to infinity in the other two directions, is removed from the rest of the field. Parallel 


_SHADOW_ SCREE! 


TURBULENCE , 
—>—t- — — ~ = $Y 











Figure 5 — Optical arrangement for taking shadow pictures. 


light after passing through this slab of the field is incident on a photographic plate (Fig. 
5). The light arriving at the plate will be more or less intense according to the distortion 
produced by density fluctuations acting as concave or convex lenses. The shadowgraph 
depends on the position of the slab relative to x; axis. A single shadowgraph gives the 
mapped field for a fixed value of z; and the complete mapping consists of a continuous 
set of shadowgraphs, one for each value of 2; . 

For small fluctuations of density, the refractive index yu is a linear function of p; 


uw = 1+ C p where C is a constant. If the photographic plate is placed closer than the 
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first focus point there are no singular points on the photographic plate, and the light 
intensity B(x, , x, 23) falling on the plate compared to the intensity, B, , of the inci- 
dent light is (7) 
Bo — B(x, , %2 , 2s) i (Ze s) 

h(z, , %2 , £3) = > TC ——- —5] G2 . 

i " : ta) B(x, ] Xo ? 3) “za-l OX; + OX2 ‘ 
We have omitted the dimensional constant in the above relation. We may mention that 
3°/ax; + 0°/dx; is the Laplacian in two dimensions and we shall later make use of this 
fact. The density field is mapped into A(x) by a simple integro-differential operator 
K(s) which of course belongs to the general class of operators considered earlier. This 
becomes clear when we express it in terms of the Dirac function and its derivatives. Thus, 


Ma) = / K(s)p(x + s) dV(s) 


where 
yo'"(s) 5(s.) + 6’’(s,) 6(s>) for ls, | < J, 
K(s) = ) 
LO for [ss | > J, 
and 6” is the second derivative of the Dirac function. The equivalence of the two ex- 
pressions for h(x) can be demonstrated by substituting for A(s) and integrating by 
parts. We now introduce the usual terminology. Let 7(&) and P(k) denote correlation 
and spectrum respectively of the density field. 


T(E) = (p(x) p(x + &))aw. 
P(k) = ws [ 7@) exp [ae-4] aV@ 


and 
Te) = [| PQe) exp [-ak-t] dV® 
P(k) is the spectral density of (p’(x)),,. in the wave number space. “Three-dimensional” 


spectrum of density, G(k), is the spectral density of (p(x)),,. with respect to wave 


number magnitude, k, 
Gk) = [ Pte) do(t) 


T(r,0,0) is the correlation usually measured by hot-wire technique. The ‘one-dimen- 
sional” spectrum, G,(k,), is the Fourier transform of 7'(r,0,0) 


rr, 0,0) = |[[ Pde) exp [iki] dk, diy dk 


Il 


/ G,(k,) exp [—tk,r] dk, 
(4.1) 


where 


G\(k,) = I P(k) dky dk, 
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If the random temperature field is statistically homogeneous and isotropic, then 


2 a 2. (f a 
Gi(ha) = 7 |] a" dls di 
(4.2) 
_1f G® », 
a an k dk. 
The inverse relation is 
G0) | 
1k) = —24| 2Gskd 4: 
G(k) 2i) a (4.3) 


Consider the mapped field A(x), since we have assumed a statistically homogeneous 
field, the statistical properties of h(x, , 22, x3) with x, and z, as variables are independent 
of z3 . Let 


BE, , 2 , O) = (A(ai , Le , Xs)h(a, + &e , Te + bo + sev. - 


This correlation can be obtained from the shadowgraph. For this purpose two slides are 
made from the shadowgraph and placed back to face. The combined pattern is the 
same as for a single slide if the position of two plates is matched. If ¢(z, , zz , 23) is the 
transparency of either plate, the transparency of two plates displaced from the matched 
position by amounts &, and &, is 

t, = (x, , Xo, Ls)t(x, + £y , Le + bo, 4s). 


Let d(x, , 2, 2%3) = to + At(z, , 22, 23) where ¢, is the average transparency. The fluctuat- 
ing transparency Ai(a, , 22, 23) = ¢,h(x, , L2 , 23) where c, is a constant, if the photo- 


rraphic process is linear. Making use of this relation and taking averages of ¢, we find 
£ I e £ 


toaw, = bo + €;(h(2; , Ye , Xa)N( Tt, + Es » Le H+ be ; La))ov. 
(4.6) 
ts + ciBlé: , & , 0). 
We can determine 6(£, , & , 0) by measuring combined transparency of two plates. 
¥() is the “auto-correlation” of A(s) 
ve) = | dss |] {8’%(s,) 4(s,) + 8's) 4(s,)} 
,6’"(s, “+ 7;) 6(S» + To.) + 5’'(8. + T>) (s, > ie 7,)} ds, ds, (4.7) 
= {8'"(r,) 8(7.) + 6'"(7.) 6(7,;) + 28’"(7;) 6/72) }(2l — | 73 |) | 73 | < 2. 
S(k) is the power sensitivity spectrum of K(s) 
S(k) = ¥(<) exp [—7k-] dV(7) 
(4.8) 


= 403 + ay (bly 


“3 


The shadowgraph method responds to second derivatives of the density field, this 


accounts for the factor (k? + 3)”. Since the light passes through a path 21 (which 
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corresponds to a hot-wire of length 2/ in the previous example of velocity spectrum 
measurement) the factor 4 [(sin k,l)/k3]’ appears in the sensitivity spectrum. If 


1 . 
rie) = 4s [ ae) exp [ak-4] av) (4.9 
then 
T(k) = S(k)P(k) 
(see (2.3a) 
= 4k + y(n tet iy 
ks 
for isotropic random density field 
T(k) = 4(k? + &3) (sin Bt —— ] = (4.10) 
Ks 4rk 
and 
ae) = + [ SP as + 03 (sin ta! ks uty exp [—ik-t] dV). (4.11) 
J ks; 


If 7 is much larger than any scale of turbulence, then [(sin k,l)/k;]° acts almost like a 
Dirac function, i.e., we can replace G(k)/k? by {G[(kKi + k3)4]}/(k? + k) and integrate 
with respect to k; . Under this assumption 


Pb 9) = ff clas + KDE + HD) exp [aks + kebd)] dey dy. (4.12) 


Introducing polar coordinates 


&, = Tr cos@¢; k, =v cos y: 

f, =rsing; kz =v sin y. 
We have 

Bir) = 2nl | v°G@)Solvr) dv, (4.12a) 

/o 

and the inverse relation 

GO) = 5 | 7800) Solr) ar. (4.12b) 

2rl /0 


The measured correlation 6(r) and v’ G(v) are Fourier-Bessel transforms of each other 
(8). For » << 1» G(v) ~ »* making the first and the third moments of 8(r) always zero. 

The relation between the conventional correlations 7'(&) and the correlation of the 
shadowgraph is given by equation (2.2) 


a) = [ W(2)T& - 2) aV(®). 


Substitute ¥(*) from equation (4.7) and integrate with respect to'r; and 7, 


B,(&) = [. (21 — | 73 (S42 +2 ae + 2 )re 73) ATs. 
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For isotropic density fluctuations 


2 + 15)"] drs, 


ctr 


Bg + 8) ] = 2] Ql — r)ATIE + 
“0 


| (4.13) 


2] (21 — 1)A27(@? + 73)?] drs, 


*@Q 


B(r) 
where A: is the bilaplacian in two dimensions for circularly symmetric case. If 1 is much 
larger than any scale of density fluctuations, then 

B(r) ‘eid 271 (2 2\1/2 é 
7 =4 | ART G? + 73)? ] drs , (4.13a) 


= A’g(r). (4.13b) 


r’ log r is the elementary solution of bilaplacian in two dimensions (9). If we assume 
that r(r),.. ~ 1/r‘; « > 0 then the solution of (4.13b) is 


=) a2 


d(r) = = | | EB(E)(r? + 2? — 2ré cos 6) log (r? + & — 2ré cos 0)'” dédé 


a 


ai | BELG? + &) log r + £7] dt + +; | BOE’ + £) log — + 1°] dé. 


We note that the first and third moments of 8(£) are zero, using this fact: 


or) = + | a@elo? + 8) logr + F]de+ 3 | Belo? + #) loge + 1°] dg 


“r 


substitute for ¢(r) from 4.13a 


g¢(r) = 4 | T[(r° a 73) "| dr. 





(4.13¢) 
= 7 | B(e)El(r? + #)(log E — log r) + r? — &"] dé. 
This integral equation can be solved for 7'(r) (see Appendix) 
ome) a 1 [ fare) — ah(ere] <2 8 
rT(r) = On | \P(s) g (s)s] [s? ase ri? 
mt [ . ds [ B(E)E[(38° + £)(log E — log s) + 2s? — 2E"] 1 
Srl J, P si nea [s* — r*]’"” 
integrate first with respect to S. 
mg yt y ré , 7 s ds 
rT) = = | awed | (8s + Flog § — log 8) + 28° — 2] eam, 
(4.14) 


— | pene, r) de, 


T(r) Srl . 


II 
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where 
Mg,7) = —5 [Be + Plog e — log s) + 2s — 22) 2S, 
_ * s =f) 


(4.15) 


|—-3V7P —r. 


2” + 2 -1 
= E ae cos 


From (4.13c) and (4.14) we see that scale of density fluctuations Z, and the micro- 
scale \, can be expressed directly in terms of B(£). 


L, = | T(r) ar / TO = “2 | B(é)E" log é at / | B(E)E* dé, 


vr is 


and 


at te -2r(0) / 70) = —6 | BEE az / | 


We denote by y,(k,) the “one-dimensional” spectrum of the mapped field. 


B(E) dé. 
v,(k,) = I| I'(k) dk, dk; 


= / Pe? 4 ye (sin fal GO) on. dls 
== || (ki + a dks dks . 
For the case of “infinitely long” light path, we can simplify this expression (see (4.11) 


and (4.12)) 


Gu > | (ki + k2G(K + 3)'?] dk 
(4.16) 
~f row 
~ I WG) Gay 


The “three-dimensional” spectrum can be recovered from the measured spectrum by 


solving equation 4.16 (see Appendix). 


WE). a 2 fey —K oe __ - 
Gen. = Fie [ (nO) Garage (4.17) 
Also 
: es 
G(k,) = i G(k) 
(4.18) 
_ 1 fd f° _lev@Ve 
7 ml I I? I (¢° _ k’)' . dé, 
integrate first with respect to k 
- —1 7° ff m#-# 3, §+¢- x)" 
Oth) «2 Aa ee ~ os he eC 
=>] nae py — qe lon 3 dt. (4.19) 


We can also use Schlieren instead of shadow method to measure correlation and spectrum 
of random density field. If we assume that the light path is larger than scale of density 
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fluctuations then true correlation and spectrum of density fluctuations can be recovered 
from the measured quantities in the same way as in the case of shadowgraph. 

5. Concluding remarks. The above considerations are useful in many other prob- 
lems besides the problem of measurement of random fields. 

In the statistical treatment of a random field it is sometimes useful to formulate the 
problem in another field which is derived from the original field in a specific manner. 
The reason for this may be theoretical convenience or that the mapped field is much 
more suitable for making physical hypotheses. 

The infinite hot-wire or light path corresponds to the case where we have averaged 
a three-dimensional random field in one direction (along the wire length) and the result 
is a two-dimensional random field. We have shown that no information is lost in the 
homogeneous case if the dependence of correlations on the distance along the wire length 
can be expressed in terms of the vector perpendicular to this direction. Having clarified 
the case of the infinite light path, we see that the shadowgraph method is essentially 
the determination of the field due to random “charge” distribution, i.e., we know the 
spectrum and correlation of the charge distribution and we want to find the spectrum 
and correlation of the “potential’’. 

In the case of a dissipative medium we can look at the field for positive times as the 
mapping of the field at time ¢ = 0 by a time dependent operator. 
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Appendix 
Solution of the integral equations 
>. [ —_o&) de ' 
va) = | ae ay (Al) 
p(t) = | g{(x? + &)'"] dé. (A2) 


<0 


Equation (A2) can be reduced to (Al) by the transformation y = (x? + ¢”)!. For the 
solution of (Al), let — = z ese 6 then 


1 ax/2 
¥(x) = - | o(x esc 6) dé. 
This equation can be reduced to Schlémlich equation by the transformation x = 1/y 


let 


v/y) _ ) (1) = 
“pile Vv(y) and = @ “Sa P(y), 


aw/2 


V(y) = | &(y sin 6) dé. 
Solution of Schlémlich equation is (10) 


w/2 
ay) =2 WWlo+2y[ Wysin 6 ae, 


T T /0 
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or in terms of the original variables 


2 :; 2s f”" 2 
d(x) = | — ry¥(4) a [ [x ese OY(x ese @)]’ ese 6 dé 
T r=@ T 


“0 





(A3) 
Of gare Ed—= 
_ e / [Ev(E)] (¢? na x’)? 
Alternatively, 
,. [| —o@ dé 
¥(x) = | Ee (sre 2)? 
can be solved by the transformation 
é — (é,) wi” 
z = (z,)"”, 
o(&) = g(t.) '”"] = &€), 
W(x) = v[(x,)~"*”?] = V(z,). 
Substituting these in our equation we get Abel integral equation. 
Wa) 1p? BG) a 
(z,)'”" 2 /0 [é, (x, en §,)]'"" 
Its solution is (10) 
20(t) __1 d [" WE) dé, _ 
(z,)"" rdz, Jo [E, (x, = i? 
or in terms of original variables 
wv. ft ie... 
o(z) = x dx ah ee — a*)'”"" (Ad 
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—NOTES— 
PURE BENDING AND TWISTING OF THIN SKEWED PLATES* 


By ERIC REISSNER (Massachusetts Institute of Technology) 


1. Introduction. The following is concerned with the problems of pure bending and 
twisting in the theory of transverse bending of thin plates. Known results for rectangular 
plates of uniform thickness will be extended to skewed plates. 

It was shown by Kelvin and Tait that the problem of St. Venant torsion of a thin 
rectangular plate with edges parallel to axes x and y has the solution w = @ry where w 
is the deflection function of the plate and @ the (constant) angle of twist per unit length. 
Within this theory, which neglects transverse shear deformation, the torque is applied 
to the plate by means of concentrated forces of suitable direction which act at the corners 
of the plate. 

We will show here that Kelvin’s and Tait’s solution is readily extended to skewed 
plates of uniform thickness. In so doing we obtain, in particular, the influence of the 
angle of skew on the torque-twist relation for the plate. We also obtain a solution for 
the problem of pure bending of a skewed uniform plate. We find that pure bending of 
the skewed plate is associated with a twisting deformation the relative magnitude of 
which depends on the angle of skew. 

2. Formulation of the problem. Let x, y be mutually perpendicular directions in the 
undeflected middle surface of the plate. The differential equation for a uniform plate bent 
by edge forces and moments only is of the form 

V’°V*w = 0 (1) 
where w is the deflection of the plate. Stress couples MW, , WM, and M,, , defined in the usual 


way, are 


M, = —D(é’w/dx* + vd’w/dy’), 
M, = —D(ew/dy? + vd*w/dz’), (2) 
M,, = —(1 — »—)D &w/dx dy. 


For what follows we have no need for the corresponding expressions for the transverse 
stress resultants. 
We consider plates bounded by two straight lines z = +/ and by two straight lines 
y +4e — x tan A. The angle A is the angle of skew of the plate, 2/1 is the span of 
the plate and c is the chord of the plate. 
Expressions for bending moment /, and twisting moment ./,, acting along the edges 
= +c — x tan A follow by means of the usual transformation formulas for plate 


Y i> 
bending couples, 
M, = M, cos’? A + M, sin’ A + 2M,, cos Asin A, 
(3) 
M,, = (M, — M,) cos Asin A + M,,(cos’ A — sin’ A). 


Received February 18, 1952. 
*The present paper is a report on work done under the sponsorship of the Office of Naval Research 
under Contract N5-ori-07834 with Massachusetts Institute of Technology. 
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In addition to this we need Kelvin’s and Tait’s result that there occur concentrated 
forces P at the corners of the plate given by 


P = (M,, + |? (4) 


and that the effective transverse edge stress resultant FR, is of the form V, + 0M,,,/ds. 
3. Choice of deflection function. We shall find that for the problems of twisting 
and bending which are here considered it is sufficiently general to assume a deflection 


function 
w = Ax” + 4By’ + Cry (5) 
where A, B and C are constants. We then have a uniform distribution of stress couples 
M, = —D(A + vB), M, = —D(B + vA), M,, = —(1 — »)DC, (6) 


vanishing transverse stress resultants V, and V, , and vanishing effective edge stress 


resultants R, . 
4. Twisting of skewed plates. The following boundary conditions must be satis‘ed 


x= +l; M, = 0; (7) 
y = +j$c — xtan A; M, = 0. (8) 

In addition to this we have that the applied torque T is given by 
T = cP. (9) 


Equations (7), (8) and (9) are three simultaneous equations for the determination of 
the three constants A, B and C. We obtain 


‘ A T tan A 
—— : B= T tan i eee. (10) 
2(1 — »)cD (l—-s ed’ (1 — »)cD 
and therewith 
(y E tan A 2 9 | 
y= "aed =p). i. (11 
o= “ST aD ise” ~*) ) 
The meaning of this result becomes somewhat more transparent if we introduce a 
new chordwise coordinate 7 counted from the centerline y = —z tan A of the plate, by 
setting 
y =n — xztan A. (12) 
We then have 
T 2 tan® A tanA . (1+ tan” A)tanA , 4 
w = (1 a eg — SS -S ma}. (13) 
~ 21 — vceD 1+» l+p 1+» 
We may define an effective angle of twist 6 per unit length by means of the expression 
§= (w), e/2 (W)ya—e/2 (14) 
cx 


and (14) gives the following result for this effective angle of twist 


Combination of (13) « 

T ( 2 tan” 4) ip 

=-— —— {1 + ——— ’ li 
. 2(1 — »)cD > v 1 
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Equation (15) indicates that the skewed plate has a smaller torsional rigidity than 
the unskewed plate and in which way the torsional rigidity varies with the angle of 
skew <A. 

5. Pure bending of skewed plate. We denote the applied moment by M. The 
following boundary conditions must be satisfied 


z=+l, cM, = M, (16) 

y = +3c — xtan A; M, = 0. (17) 

In addition to this we have the condition of vanishing torque, or of vanishing corner 
forces P, which becomes 

M,, + M,, = 0. (18) 


From equations (16) to (18) we obtain the following expressions for the coefficients A, 
B and C in w 


M vM ‘ M tan A (19) 





(1 — cD’ hae Y*)cD’ ‘21 — veD’ 
The deflection w is now 


M 


, = —_ — — 2 — € 
u a — cD [x? — vy (1 + ») tan A zy]. (20) 
In terms of the chordwise variable n defined by (12) this becomes 
ore _ — J vm , eo 2 
ea PeD {(1 + tan’ A)z (1 — ») tan A an — v7 j. (21) 


We see from (21) that the skewed plate has a smaller bending stiffness than the 
unskewed plate and that moreover the applied bending moment produces also a torsional 
deformation. 


ON THE ERROR TERM IN INTERPOLATION FORMULAS 
By IVAN NIVEN (University of Oregon) 


From n + 1 known values of a function f(z) an approximating polynomial p(x) 
of degree n can be obtained. If f(x) and its derivatives to order n + 1 are continuous, 
then the error term can be calculated from 


inne | I 
f(z) — pa) = in + 1)! (x — 2), 
where f(z;) are the known values of the function for j = 0, 1, --- ,n,z is an arbitrary point 
in the interval (x) , z,) with z # 2; for every j, and & is some point in the interval de- 
pendent on z. The purpose of this note is to indicate that no limit on the error can be 
obtained by replacing the term f*” (é) by some function of an n-th difference A”**y. 
(This replacement is made, for example, in Numerical Mathematical Analysis by J. B. 











398 NOTES [Vol. X, No. 4 
Scarborough, 2nd edition, pp. 99-103, the mistake arising from identification of two 
é-values in the interval which may be distinct.) 

Consider the given data f(j) = j for 7 = 0, 1, --- , n, for which the polynomial 
approximation is p(x) = x. Can the error be bounded by any function of 2 , --- , 2, , 
Yo = f (2), , Yn = f(x,)? That it cannot is clear from the function f(z) = «+k 
sin rz. We have f(3) = 4 + k, p(4) = 3, and f(3) — p(4) = k. Since k is independent 
of z; and f(z,), this error value cannot be bounded by any function of these without 


strong hypotheses on the function f 


A NOTE ON MY PAPER 


ON AXIALLY SYMMETRIC FLOW AND THE METHOD 
OF GENERALIZED ELECTROSTATICS* 


QUARTERLY OF APPLIED MATHEMATICS, 10, 197-213 (1952) 


By L. E. PAYNE (University of Maryland) 


It has been brought to the author’s attention that the flow problem for a spindle was 
considered by E. W. Hobson (‘‘On a class of spherical harmonics of complex degree with 
application to physical problems’’, Trans. Camb. Phil. Soc., 14, 211-236 (1889)). Hobson 
used a method which is entirely different from that employed by the author, but unfor- 
tunately his solution is in error. The solution to the problem is given in a corrected and 
simplified form in this paper. The only reference to Hobson’s solution which the author 
has found is given in the appendix of A. B. Basset’s ‘‘Treatise on hydrodynamics”, 
vol. 2 (1888). Basset, however, did not discuss the problem and consequently did not 


recognize the errors in Hobson’s solution. 








*Received July 21, 1952 


BOOK REVIEWS 


Finite deformation of an elastic solid. By Francis D. Murnaghan. John Wiley & Sons, 
Ine., New York, 1951. viii + 140 pp. $4.00. 


The book contains an unusually lucid exposition of the theory of finite elastic deformations, a field 
in which the basic physical ideas are often lost in the tangle of complex mathematical notations. The 
elegance of the present treatment is achieved by the consistent use of matrices. The first chapter serves 
as an introduction to vectors and matrices. The second chapter is concerned with the strain matrix, 
its behavior under transformations of the initial and fina] reference frames, its invariants, and the 
compatibility relations. The stress matrix and the general relations between stress and strain are dis- 
cussed in Chapter 4, and for non-isotropic materials in Chapter 5. Chapters 6 and 7 are devoted to the 
application of the theory to specific problems such as simple shear, simple tension, or torsion of a circular 
cylinder. 

Bringing clarity into a field where confusion is the rule, the volume will doubtless be welcomed by 
all students of mechanics of continua. In two respects, however, the book disappointed this reviewer. 
First, there is no reference whatsoever to the considerable volume of classical and recent work in this 
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field. While one can readily understand the author’s refusal to work his way through the maze of con- 
fusing notations and often conflicting results of these papers, it cannot be denied that a critical digest 
of at least the more important of these papers would have greatly increased the value of the present 
book. Secondly, while the basic assumptions of the theory are clearly stated, the author rarely discusses 


their physical justification. 
W. PRAGER 


Jacobian elliptic function tables. By L. M. Milne-Thomson. Dover Publications, Inc., 
1950. xi + 132 pp. $2.45. 


The book opens with 38 pages of definitions, identities, and generally useful information concerning 
the elliptic functions sn(u, m), en(u, m), dn(u, m), z(u). The first three of these are tabulated to five 
places for values of u running from 0 to 3 in increments of .01, and for values of m from 0 to 1 in incre- 
ments of .1. The function z(u) is tabulated to seven places for u ranging from 0 to 3 in increments of 


.01. The complete elliptic functions are also tabulated and en(u, 1), dn(u, 1) are given for 3 < u < 4. 


G. F. CARRIER 


Tables of the error function and of its first twenty derivatives. By The Staff of the Computa- 
tion Laboratory. Harvard University Press, Cambridge, Mass., 1952. xxvii + 276 pp. 


$8.00. 


The function (27)~!/2 {, exp (—u?/2) du and its first twenty-one derivatives are tabulated. The 
argument increments are .004 for the function and its first nine derivatives, and are .002 for the remaining 
derivatives. The argument range is from zero in all cases to 6.468, 8.236, 9.610, and 10.902. These ranges 
refer respectively to: the error function and its first five derivatives, the next five derivatives, the twelfth 
to sixteenth derivatives, and the rest of the twenty-one derivatives. 

G. F. CARRIER 


Tables relating to Mathieu functions. Prepared by The Computation Laboratory of the 
National Applied Mathematics Laboratories, National Bureau of Standards. Colum- 
bia University Press, New York, 1951. xlvii + 278 pp. $8.00. 


This book is concerned with the solutions of the Mathieu equation, uz: + (b — s cos? r)u = 0, which 


have the period 27, = 0. s is treated as a given parameter and b plays the role of the eigenvalue. The 
first thirty eigenvalues (15 corresponding to even solutions, 15 to odd) are tabulated for s ranging from 
0 to 100. Since these functions have the representation u = = D, cos nz (or sin nx), the values of D, 


are tabulated for various n, s, b,(s). Finally, the joining factors are tabulated over the appropriate range 
of s. All this is prefaced by that general information which renders the tables easy to use. 


G. F. CARRIER 


Activity analysis of production and allocation. Proceedings of a Conference. Edited by 
Tjalling C. Koopmans. John Wiley & Sons, Inc., New York, and Chapman and 
Hall, Ltd., London, 1951. xiv + 404 pp. $4.50. 


As proceedings of a conference, this book is almost necessarily a compilation of articles and papers 
of uneven quality and a superficially great diversity of topics. Nevertheless, in this reviewer's opinion it 
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portends a coming of age of economics mathematically, both in the sense of an applied mathematical 
economics, and in the sense of a discipline whose techniques may be employed profitably in some older 
‘‘well-mathematized”’ areas. The crux of this is the brilliant theory of interacting productive aggregates 
of T. C. Koopmans (Chapter III) which furnishes a background for economic studies comparable in 
scope and utility to the linear steady state theory of electrical circuits in its realm. 

Although the developments of this volume have profound implications for centralized planning, 
allocation, decentralization of authority, and orientation of research (profit-wise) in industry, they 
were largely stimulated by practical problems of military logistics. Thus the Army Air Forces group 
under M. K. Wood and G. B. Dantzig developed methods for (a) programming interdependent activities, 
e. g., the Berlin Airlift, (b) analysis of inter-industry commodity flows (generalizing and employing 
W. W. Leontief), (c) effective computation. 

Mathematically, the central problem—the “linear programming” problem—is maximization of a 
linear functional (sometimes a vector) of non-negative variables subject to linear inequalities, e. g. the 
“transportation problem’, to ship specified quantities of a product to n destinations from m origins 
under specified availability restrictions and travel costs so as to have least total shipment cost (resp. to 
find the “efficient”’ combinations of activities in production). Requisite theory of convex cones is devel- 
oped ab ovo by Gerstenhaber, Gale and Arrow with interesting extensions by Gale, Kuhn and Tucker in 
the direction of equivalent maximal and minima 





principles. There is close connection with 2-person game 
theory—ev 





y game can be written as a linear programming problem, and vice-versa usually, so that 
results and techniques in either bear on the other 

Computationwise (and for some theoretical questions) in linear programming, the finite step “‘sim- 
plex’? method of Dantzig, which explicitly informs of arrival at a maximal solution, is one of the most 
effective tools yet developed—far outclassing, say, the Dines elimination technique. Heuristically em- 
ployed up to the present, since Dantzig’s proofs apply only to extremely special situations, recent results 
of the reviewer (‘‘Mathematical Background of Linear Programming”, Carnegie Tech-Air Forces Project 
in Intra-Firm Analysis, July 1951 provide for its applicability to every situation. 

To supplement the 1 i 





eadingly meager specific applications of linear programming detailed, current 
ones include: inversion of matrices, solution of simultaneous linear equations, optimal award of contracts 
to bidders, blending of aviation gasolines, and even plastic limit design of structures. 


A. CHARNES 


Thermodynamics of irreversible processes. By S. R. DeGroot. North-Holland Publishing 
Company, Amsterdam, and Interscience Publishers, Inc., New York, 1951. xvi + 
242 pp. $4.00. 


} 


The subject of irreversible processes has developed in relatively complete form within the last 
decade. This book shows how the ideas of entropy flow and entropy production together with the Onsager 
reciprocal relations lead to a theory of the thermodynamics of irreversible processes. The first two 
chapte rs of the book give an account of the theory, especially the ¢ nsager relations; the following seven 
chapters include many examples from physics and chemistry and, among other topics, discussion of the 


application of the theory to heat conduction, electrical conduction (with and without a magnetic field), 


relaxation phenomena in continuous single component systems, discontinuous systems (with and without 
chemical reactions), ordinary diffusion, thermal diffusion, viscosity, diffusion potentials, thermoelec- 
tricity, reaction rates, electrochemistry, electrokinetic effects, and interference of a chemical reaction and 
a rel 





xation phenomenon. The last two chapters of the book are concerned with special topics such as the 
definition and discussion of stationary states of various order and further discussion of the foundations 
of the Onsager relations and the general theory 

The author has marked sections of the book so that it can be read in three cycles, the first for a gen- 
eral idea of the thermodynamics of irreversible processes, the second includes the examples but omits the 
statistical basis of the theory, and the third includes the entire monograph. It seems clear that the book 
will be tremendously useful to most physical chemists and to many physicists. The book assumes a 
background in thermodynamics and physical statistics although much of it can be read without the latter. 


Roun TRUELL 
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Die Welt der Vektoren. Finfithrung in Theorie und Anwendung der Vektoren, Tensoren 
und Operatoren. By Franz Ollendorff. Springer-Verlag, Wien, 1950. viii + 470 pp. 
$9.50. 


Dr. F. Ollendorff, other books by whom the reader has certainly used and enjoyed, is at present 
Professor at the Hebrew Technical College, Haifa (Israel). The present work is derived from a course of 
vector analysis, preparatory to modern physics, given to the students in their senior year. Its purpose is 
thus frankly utilitarian, but the book is penetrated through and through with the author’s admiration for 
the beauty of its subject. This combination is so effective that the reader is often ashamed of voicing his 
objections to himself, as if he were interrupting a beautiful lecture. 

The book contains eight chapters. Chapters 1, Vectors and Scalars, and 2, Vector Fields, contains the 
standard material and its application to geometry, mechanics, fluid mechanics and Maxwell electro- 
dynamics, together with a few less common examples. 

Chapters 3, Vector Analysis in A ffine Space, 4, Tensor Algebra, and 5, Tensor Analysis in A ffine Space, 
present the fundamental tools of tensor analysis applied to an affine space of z dimensions, with metric. 
(The possible absence of a metric is not discussed). There is an interesting treatment of pseudo-scalars 
ind of axial vectors of both kinds, inspired by Brillouin but really distinct, and connected with the differ- 
ent definitions of the positive side of a surface element in Stokes’ and Gauss’ theorems. Chapter 6, 
Minkowski Space, applies this material to special relativity, as well as to de Broglie and Schrsdinger 
Waves. 

Chapter 7, Riemannian Space, gives the essentials of parallel displacement, geodesics and covariant 
differentiation. This is then used in a short but beautifully clear presentation of the dynamics of general 
relativity. 

Chapter 8, Hilbert Space, generalizes from real to complex scalars and from a finite to an infinite 
number of dimensions. The mathematical reader, who will have had misgivings during the treatment of 
tensor analysis, will here become quite restive. The author however keeps his goal firmly in mind, and 
rapidly reaches the pass from which he can lead his audience into the green pastures of linear integral 
equations, matrix mechanics, Hamilton and Pauli operators. 

A feature of the book is its richness in physical examples, to which it is impossible to do justice in a 
review. As the author very well savs in his Preface, a work on vectors cannot be at the same time a text- 
book on Physics. Still every writer on Vector Analysis is anxious to give examples, and some have given 
such a condensed treatment of several extended physical fields that they end by serving neither the 
beginner nor the advanced reader. It seems to us that Prof. Ollendorff has struck on this point a very 
happy balance. The treatment of the physical examples (with the exception perhaps of some at the very 
end) is sufficiently extended to give the reader, due in part no doubt to the author’s skill, no feeling of 
undue condensation or hurry. Indeed any well-prepared student should find this book a most readable 
text, and it will introduce him to some of the most beautiful questions of mathematical physics beside 


those from his own special field. The Haifa Technical College should assuredly be proud of such a course, 


ot the le eturer and ol his book. 
P. Le CorBEILLER 


Theory of perfectly plastic solids. By William Prager and Philip G. Hodge, Jr. John Wiley 
& Sons, Inc., New York and Chapman « Hall, Limited, London, 1951. x + 264 pp. 


$5.50. 


The student and research worker of today is in a fortunate position, as compared with his predecessor 
a generation ago, when he wants to penetrate a new field of science without turning to the scattered 
original literature. Instead of comprehensive text books, covering whole subjects such as Hydrodynamics 
or the Theory of Elasticity, he may choose among a great variety of more specialized monographs which 
make it possible for him to avoid devouring a lot of unnecessary material. The accompanying disadvan- 
tages of this process of development are obvious, but are probably overestimated by the older generation 
of scientists. The development has been actively promoted by American publishers. 

The book under review forms an outstanding example of this type of modern monographs. It pre- 
tends to be a first introductory treatment, written on an intermediate level, and should be easily under- 
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stood by the better students, having taken the basic courses in mathematics and mechanics of most 
engineering schools. With this principle in view, it might be said that the complementary mathematica] 
treatment, now given in appendices, as well could have been incorporated in the main text, at least from 
the point of view of a European reader. To facilitate the penetration of the concentrated subject, ample 
references to elementary textbooks, even with indication of page numbers, are given. 

Faithful to the scope covered by its title, the book gives a rather complete treatment of the theory of 
perfectly plastic solids according to v. Mises or Prandtl-Reuss. Interesting applications to engineering 
problems are given. The chapters are: 1. Basic concepts, 2. Trusses and beams, 3. Torsion of cylindrical 
or prismatic bars, 4. Plane strain: problems with axial symmetry, 5. Plane strain: general theory, 6. Plane 
strain: specific problems, 7. Plane strain: contained plastic deformation. Limit analysis, 8. Extremum 
principles. Cartesian tensors with subscript notation and summation conventions do appear only in 
chapter 8. Each chapter contains exercises and a list of references. The text and figures are distinguished 
by their clarity and stringency. 

Those who look for physical] aspects of the theory of plasticity, such as details as to recent develop- 
ment of the flow theory as compared with the so-called deformation theory of plastic solids, are referred 


to other publications. 
Fo.kE K. G. Opqvisr 


Operational calculus based on the two-sided Laplace integral. By Balth. Van Der Pol and 
H. Bremmer. University Press, Cambridge, 1950. xiii + 415 pp. $10.00. 


The material in this book evolved from a series of lectures by the authors in the period from 1938 to 
1940. Written up and extended during the German occupation of the Netherlands, the text was later 
translated from Dutch to English. 

The book is concerned with the application of operational calculus to mathematics, physics and 
engineering, and the stated aim is a modern treatment of the subject. 

The chapter headings are as follows: I. General Introduction, II. The Fourier Integral as Basis of 
the Operational Calculus, III. Elementary Operational Images (transforms), IV. Elementary Rules, 
V. The Delta of Impulse Function, VI. Questions Concerning the Convergence of the Definition Integral, 
VII. Asymptotic Relations and Operational Transposition of Series, VIII. Linear Differential Equations 
with Constant Coefficients, IX. Simultaneous Linear Differential Equations with Constant Coefficients; 
Electric Circuit Theory, X. Linear Differential Equations with Variable Coefficients, XI. Operational 
Rules of More Complicated Character, XII. Step Functions and Other Discontinuous Functions, XIII. 
Difference Equations, XIV. Integral Equations, XV. Partial Differential Equations in the Operational 
Calculus of One Variable, XVI. Simultaneous Operational Calculus, XVII. Grammar (Operational 
Rules), XVIII. Dictionary (of transforms). 

There are a number of excellent features of this book. One chapter, and other sections of the book, 
on the treatment of the delta function seem particularly good; they deal at some length with various 
features of the delta function, including applications in connection with differential and integra] equa- 
tions. The impulse function is formulated in terms of the Stieltjes integral; the transform of the delta 
function is discussed; functions approximating the delta function are considered. The use of the impulse 
function in obtaining the Green function of inhomogeneous differential equations is discussed. The sections 
treating integral equations, partial differentia] equations, and the transformation of functions of more 
than one variable by means of multiple Laplace integrals are also interesting, for one reason, because of 
the use of the delta function techniques in connection with these topics. 

The treatment of the operational calculus given in this book is based on the two-sided Laplace trans- 
form, i.e., with limits — ~ and = instead of the usual limits 0 and ~. This treatment demands that the 
separate one-sided integrals must have a common region of convergence otherwise the two-sided Laplace 
transform cannot exist. The overlapping of the convergence regions or the strip of convergence must be 
ascertained and may require somewhat more care and attention than the usual one-sided integrals. 
Explicit formulae for the convergence strip for any given original function h(t) are obtained and discussed 
in the chapter on convergence of the definition integral. 

The book seems to be very carefully written and it should be of real value to the student of operational 


methods. 
Roun TRUELL 























